This code generates the supplemental figures from the outputs of the model.
binder, where is the .yml file for the conda virtual environment.
data, where are the outputs of the model NEMO-LIM3, the outputs of the model MITgcm, and the observations from the Qikiqtarjuaq ice camps 2015 and 2016.
figures_supmat_progress, where are the supplemental figures.
This notebook was run in a conda virtual environment. The following commands will automatically create the environment and launch the Jupyter notebook:
cd <my_path>/wintertime/
conda env create -f binder/environment.yml
conda activate wintertime
cd utils/python/MITgcmutils
python setup.py bdist_egg
cd build
pip install MITgcmutils
cd ../../../..
jupyter notebook
Once the environment wintertime is created, the Jupyter notebook can be launched simply with:
cd <my_path>/wintertime/
conda activate wintertime
jupyter notebook
import datetime
import gsw
import gudinfo
import matplotlib as mpl
import matplotlib.cm
import matplotlib.colors
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import netcdf_tools
import numpy as np
import pandas as pd
import plotly.graph_objects as go
import plotly.express as px
import read_mitgcm
import scipy
import scipy.interpolate
import statsmodels.formula.api as smf
import statsmodels.nonparametric.smoothers_lowess
import TaylorDiagram
import vstats
plt.close("all")
years=range(1,11)
nbyears=len(years)
first_year=np.arange(0,365)
last_year=np.arange(365*(nbyears-1)+0,365*(nbyears-1)+365)
array1d_iT1y_iT=last_year
first_year366=np.arange(0,366) # first year for heatmaps
drF is the r cell face separation, meaning the thickness of each depth layer (in m)
It corresponds to delR on https://mitgcm.readthedocs.io/en/latest/getting_started/getting_started.html#grid
gridfile='data/DataS8_output_mitgcm/exp0/grid.t001.nc';
drF=netcdf_tools.read_netcdf(gridfile,'drF')
RC is the r coordiante of cell center (in m). We switch sign to have positive depths.
RC=-netcdf_tools.read_netcdf(gridfile,'RC')
RF is the r coordinate of cell interface (in m). We switch sign to have positive depths.
RF=-netcdf_tools.read_netcdf(gridfile,'RF')
RF_above151=np.array(RF[RF<151])
See a description of the sea ice camps in Mundy et al. (2014) in Mar. Ecol. Prog. Ser. DOI: https://doi.org/10.3354/meps10587.
The file Ardyna_Chl_Resolute.csv contains $in\ situ$ chlorophyll $a$ at the Resolute sea ice camps 2010 and 2011 ($74.71^\circ N$, $-92.25^\circ E$). Variable 'chl_mg_m2' is the vertically integrated (0-100m) chlorophyll $a$ measured by fluorimetry (Ardyna et al., 2020, in $Elem.\ Sci.\ Anth.$ https://doi.org/10.1525/elementa.430). We acknowledge CJ Mundy for the collection of data.
def load_chl():
fname='obs/Ardyna_Chl_Resolute.csv'
df = (
pd.read_csv(
fname,
dtype={
'Date': str,
'chl_mg_m3': np.float32,
'chl_mg_m2': np.float32,
},
parse_dates=['Date']
)
.rename(
columns=dict(
chl_mg_m2='chl'
)
)
)
df['doy']=df['Date'].apply(lambda x:x.timetuple().tm_yday)
return df
chl_df=load_chl()
chl2010_df=chl_df[chl_df.Date.dt.year==2010]
chl2011_df=chl_df[chl_df.Date.dt.year==2011]
in $mg\ Chla\ m^{-2}$
nT=365
array1d_iT_chl2010=np.full(nT,np.NaN)
for doy in chl2010_df['doy']:
iT=doy-1
chl=chl2010_df.loc[chl2010_df.doy==doy,'chl']
array1d_iT_chl2010[iT]=chl
nT=365
array1d_iT_chl2011=np.full(nT,np.NaN)
for doy in chl2011_df['doy']:
iT=doy-1
chl=chl2011_df.loc[chl2011_df.doy==doy,'chl']
array1d_iT_chl2011[iT]=chl
blooms_df=pd.DataFrame(index=[2010,2011],
columns=['tI','cI','tP','cP','m'],
dtype=np.float32)
The bloom initiation was identified when the accumulation rates ($d^{-1}$) were positive for one week (Boss and Behrenfeld, 2010):
accumulation_rate=$\frac{chl_{iT}-chl_{iT-1}}{chl_{iT-1}}$
$In\ situ$ data was smoothed with LOWESS (locally weighted scatterplot smoothing) and interpolated before calculating the accumulation rates.
x=first_year
y=array1d_iT_chl2010
df=pd.DataFrame({'Day of year':x,'Observed Chl a vert. integrated (0-100m)':y})
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,700])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
Figure. Scatterplot of day of year and observed Chl a vert. integrated (0-100m).
y_hat=statsmodels.nonparametric.smoothers_lowess.lowess(y,x,frac=1/4)
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Add the prediction line
fig.add_traces(go.Scatter(x=y_hat[:,0], y=y_hat[:,1], name='LOWESS, frac=1/4', line=dict(color='red')))
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,700])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
# ------- Define interploation functions -------
# Linear - draws a line between the two nearest points and calculates y value based on the slope of that line
f_linear = scipy.interpolate.interp1d(y_hat[:,0],
y=y_hat[:,1],
bounds_error=False)
# Create a new set of points with x values
indices_not_nan=np.argwhere(~np.isnan(y)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
xnew=np.empty(end-start+1)
for iT in range(start,end+1):
xnew[iT-start]=iT
# Find y values based on the interpolation method
array1d_iT_obsvintchlHPLClowess=f_linear(xnew)
start_edge=np.empty(start)
start_edge[:]=np.NaN
end_edge=np.empty(365-1-end)
end_edge[:]=np.NaN
array1d_iT_obsvintchlHPLClowess=np.concatenate((start_edge,
array1d_iT_obsvintchlHPLClowess,
end_edge))
x=first_year
y=array1d_iT_obsvintchlHPLClowess
df=pd.DataFrame({'Day of year':x,'Smoothed and interpolated observed Chl a vert. integrated (0-100m)':y})
indices_not_nan=np.argwhere(~np.isnan(array1d_iT_obsvintchlHPLClowess)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
# no NaN between first and last smoothed observed value
assert np.array_equal(indices_not_nan,np.arange(start,end+1)), (
'There is at least one NaN in array1d_iT_obsvintchlHPLCsmooth before '
'its first and last valid values.')
array1d_iT_r2010=np.empty(365)
array1d_iT_r2010[:]=np.NaN
for iT in range(start+1,end):
chliT=array1d_iT_obsvintchlHPLClowess[iT]
chljT=array1d_iT_obsvintchlHPLClowess[iT-1]
acc=(chliT-chljT)/chljT
array1d_iT_r2010[iT]=acc
x=first_year
y=array1d_iT_r2010
df=pd.DataFrame({'Day of year':x,'Accumulation rates':y})
found=False
iT=0
while not found and iT<=364:
array1d_iT_accweekbefore=array1d_iT_r2010[np.arange(iT-6,iT+1)]
if(np.all(~np.isnan(array1d_iT_accweekbefore)) \
and \
np.all(array1d_iT_accweekbefore>0)):
found=True
else:
iT=iT+1
tI=iT
year=2010
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tIs=datetime.datetime(year,1,1) + datetime.timedelta(tI)
Day of year
tI=tI+1
tI
140
tIs.strftime("%B %d")
'May 20'
cI=array1d_iT_obsvintchlHPLClowess[tI-1]
(in $mg\ Chl\ m^{-2}$)
cI
10.521742105483984
blooms_df.loc[2010,['tI']]=tI
blooms_df.loc[2010,['cI']]=cI
with plt.style.context('mplstyles/exp.mplstyle'):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
fig,axs=plt.subplots(2,1,sharex=False,figsize=(16, 6))
# --- ACCUMULATION RATE
ax=axs[0]
h1=ax.plot(first_year,array1d_iT_r2010,'-',color='blue',lw=3)
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks(np.arange(-0.5,1.0,0.5))
ax.set_ylim(-0.5,0.5)
ax.set_ylabel('Accumulation rate ($\mathrm{ d^{-1} }$)')
ax.plot(xlims,[0,0],'k--',color='gray')
# adding the vertical line at start of bloom development
ax.axvline(tI-1,color='green',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[1]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.12,0.15 ,0.65,.60])
axs[1].set_position( [0.12,0.10 ,0.65,.01])
x=first_year
y=array1d_iT_chl2011
df=pd.DataFrame({'Day of year':x,'Observed Chl a vert. integrated (0-100m)':y})
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,700])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
Figure. Scatterplot of day of year and observed Chl a vert. integrated (0-100m).
y_hat=statsmodels.nonparametric.smoothers_lowess.lowess(y,x,frac=1/4)
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Add the prediction line
fig.add_traces(go.Scatter(x=y_hat[:,0], y=y_hat[:,1], name='LOWESS, frac=1/4', line=dict(color='red')))
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,700])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
# ------- Define interploation functions -------
# Linear - draws a line between the two nearest points and calculates y value based on the slope of that line
f_linear = scipy.interpolate.interp1d(y_hat[:,0],
y=y_hat[:,1],
bounds_error=False)
# Create a new set of points with x values
indices_not_nan=np.argwhere(~np.isnan(y)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
xnew=np.empty(end-start+1)
for iT in range(start,end+1):
xnew[iT-start]=iT
# Find y values based on the interpolation method
array1d_iT_obsvintchlHPLClowess=f_linear(xnew)
start_edge=np.empty(start)
start_edge[:]=np.NaN
end_edge=np.empty(365-1-end)
end_edge[:]=np.NaN
array1d_iT_obsvintchlHPLClowess=np.concatenate((start_edge,
array1d_iT_obsvintchlHPLClowess,
end_edge))
x=first_year
y=array1d_iT_obsvintchlHPLClowess
df=pd.DataFrame({'Day of year':x,'Smoothed and interpolated observed Chl a vert. integrated (0-100m)':y})
indices_not_nan=np.argwhere(~np.isnan(array1d_iT_obsvintchlHPLClowess)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
# no NaN between first and last smoothed observed value
assert np.array_equal(indices_not_nan,np.arange(start,end+1)), (
'There is at least one NaN in array1d_iT_obsvintchlHPLCsmooth before '
'its first and last valid values.')
array1d_iT_r2011=np.empty(365)
array1d_iT_r2011[:]=np.NaN
for iT in range(start+1,end):
chliT=array1d_iT_obsvintchlHPLClowess[iT]
chljT=array1d_iT_obsvintchlHPLClowess[iT-1]
acc=(chliT-chljT)/chljT
array1d_iT_r2011[iT]=acc
x=first_year
y=array1d_iT_r2011
df=pd.DataFrame({'Day of year':x,'Accumulation rates':y})
found=False
iT=0
while not found and iT<=364:
array1d_iT_accweekbefore=array1d_iT_r2011[np.arange(iT-6,iT+1)]
if(np.all(~np.isnan(array1d_iT_accweekbefore)) \
and \
np.all(array1d_iT_accweekbefore>0)):
found=True
else:
iT=iT+1
tI=iT
year=2011
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tIs=datetime.datetime(year,1,1) + datetime.timedelta(tI)
Day of year
tI=tI+1
tI
145
tIs.strftime("%B %d")
'May 25'
cI=array1d_iT_obsvintchlHPLClowess[tI-1]
(in $mg\ Chl\ m^{-2}$)
cI
9.559321403503418
blooms_df.loc[2011,['tI']]=tI
blooms_df.loc[2011,['cI']]=cI
with plt.style.context('mplstyles/exp.mplstyle'):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
fig,axs=plt.subplots(2,1,sharex=False,figsize=(16, 6))
# --- ACCUMULATION RATE
ax=axs[0]
h1=ax.plot(first_year,array1d_iT_r2011,'-',color='blue',lw=3)
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks(np.arange(-0.5,1.0,0.5))
ax.set_ylim(-0.5,0.5)
ax.set_ylabel('Accumulation rate ($\mathrm{ d^{-1} }$)')
ax.plot(xlims,[0,0],'k--',color='gray')
# adding the vertical line at start of bloom development
ax.axvline(tI-1,color='green',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[1]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.12,0.15 ,0.65,.60])
axs[1].set_position( [0.12,0.10 ,0.65,.01])
The date of the bloom peak was the date of the maximum (vertically integrated) chlorophyll $a$.
tP=np.nanargmax(array1d_iT_chl2010)
tP=tP.item()
year=2010
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tPs=datetime.datetime(year,1,1) + datetime.timedelta(tP)
(in $mg\ Chl\ m^{-2}$)
cP=np.nanmax(array1d_iT_chl2010)
cP
659.6625366210938
Day of year
tP=tP+1
tP
172
tPs.strftime("%B %d")
'June 21'
Since the highest value was the last observation, the chlorophyll $a$ may have continued to increase after the last observation. It is impossible to know for sure if June 21 was really the bloom peak.
blooms_df.loc[2010,['tP']]=tP
blooms_df.loc[2010,['cP']]=cP
tP=np.nanargmax(array1d_iT_chl2011)
tP=tP.item()
year=2011
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tPs=datetime.datetime(year,1,1) + datetime.timedelta(tP)
(in $mg\ Chl\ m^{-2}$)
cP=np.nanmax(array1d_iT_chl2011)
cP
324.3826599121094
Day of year
tP=tP+1
tP
162
tPs.strftime("%B %d")
'June 11'
blooms_df.loc[2011,['tP']]=tP
blooms_df.loc[2011,['cP']]=cP
$m$ is the slope between the bloom initiation (I) and the bloom peak (P) in the log scale. A linear regression of the points between I and P was used.
obsvintchl2010_df=pd.DataFrame(data={'iT':first_year,
'obsvintchl':array1d_iT_chl2010})
obsvintchl2010_df=obsvintchl2010_df[obsvintchl2010_df.obsvintchl.notna()]
y=2010
tI=blooms_df.loc[y,['tI']].item()-1 # doy (1-based) -> doy (0-based)
tP=blooms_df.loc[y,['tP']].item()-1 # doy (1-based) -> doy (0-based)
Points after bloom initiation and before bloom peak
initiation2010_df=obsvintchl2010_df[(obsvintchl2010_df.iT>=tI)
&(obsvintchl2010_df.iT<=tP)]
Points before bloom initiation and after bloom peak. Those points were note used for the linear regression.
out2010_df=obsvintchl2010_df[(obsvintchl2010_df.iT<tI)
|(obsvintchl2010_df.iT>tP)]
loginitiation2010_df=pd.DataFrame(
data={'iT':initiation2010_df.iT,
'obsvintlogchl':np.log10(initiation2010_df.obsvintchl)})
olsres=smf.ols('obsvintlogchl~iT',data=loginitiation2010_df).fit()
m=olsres.params.iT
blooms_df.loc[2010,['m']]=m
To draw the line representing $m$.
olsx2010=[tI,tP]
olsy2010=olsres.predict(pd.DataFrame({'iT':olsx2010}))
olsy2010=np.power(10,olsy2010)
obsvintchl2011_df=pd.DataFrame(data={'iT':first_year,
'obsvintchl':array1d_iT_chl2011})
obsvintchl2011_df=obsvintchl2011_df[obsvintchl2011_df.obsvintchl.notna()]
y=2011
tI=blooms_df.loc[y,['tI']].item()-1 # doy (1-based) -> doy (0-based)
tP=blooms_df.loc[y,['tP']].item()-1 # doy (1-based) -> doy (0-based)
Points after bloom initiation and before bloom peak
initiation2011_df=obsvintchl2011_df[(obsvintchl2011_df.iT>=tI)
&(obsvintchl2011_df.iT<=tP)]
Points before bloom initiation and after bloom peak. Those points were note used for the linear regression.
out2011_df=obsvintchl2011_df[(obsvintchl2011_df.iT<tI)
|(obsvintchl2011_df.iT>tP)]
loginitiation2011_df=pd.DataFrame(
data={'iT':initiation2011_df.iT,
'obsvintlogchl':np.log10(initiation2011_df.obsvintchl)})
olsres=smf.ols('obsvintlogchl~iT',data=loginitiation2011_df).fit()
m=olsres.params.iT
blooms_df.loc[2011,['m']]=m
To draw the line representing $m$.
olsx2011=[tI,tP]
olsy2011=olsres.predict(pd.DataFrame({'iT':olsx2011}))
olsy2011=np.power(10,olsy2011)
The file Pigments, nutrients, Chlorophyll a and Phaeopigments (concentration) from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ chlorophyll a at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variables 'conc_mg_m3' is the chlorophyll a concentration measured by HPLC (mg Chl $\mathrm{ m^{-3} }$).
def load_pigment():
fname='obs/66417.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'date': str,
'date_time': str,
'latitude': np.float32,
'longitude': np.float32,
'sample_id': 'string',
'sample_code': 'string',
'sample_type': 'category',
'sample_source': 'category',
'snow_thickness': 'category',
'leg': 'category',
'code_operation': 'category',
'station': 'category',
'station_type': 'category',
'cast': np.float32,
'bottle': np.float32,
'target_depth_m': 'string',
'depth_m': np.float32,
'pigment': 'category',
'conc_mg_m2': np.float32,
'conc_ug_l': np.float32,
'conc_mg_m3': np.float32,
'filtered_vol_l': np.float32,
'vol_after_melt_l': np.float32,
'ice_added_vol_of_water': np.float32,
'ice_total_vol_water': np.float32,
'ice_calculated_vol_ice': np.float32,
'dilution_factor': np.float32,
'quality_control_peaks': 'category',
'nb_core': np.float32,
'nb_core_9cm_dia': np.float32,
'core_area_9cm_dia': np.float32,
'nb_core_14cm_dia': np.float32,
'core_area_14cm_dia': np.float32,
'area_m2': np.float32,
'description_of_sample': 'category',
'n': np.float32,
'observations': 'category',
'method': 'category',
'pi': 'category'
},
parse_dates=['date','date_time']
)
.rename(
columns=dict(
depth_m='depth',
conc_ug_l='chlfluo',
conc_mg_m3='chlHPLC'
)
)
)
df['depth']=-df['depth'] # negative depths
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
return df
chlHPLC2015_df=load_pigment()
We select the observations of chlorophyll $a$ concentrations at the ice camp 2015 measured by HPLC.
chlHPLC2015_df=\
chlHPLC2015_df[(chlHPLC2015_df.mission=='ice_camp_2015') \
& (chlHPLC2015_df.sample_type=='water') \
& (chlHPLC2015_df.pigment=='Total Chlorophyll a') \
& (chlHPLC2015_df.method=='HPLC')]
find negative chlorophyll $a$
chlHPLC2015_df['chlHPLC'][ chlHPLC2015_df['chlHPLC']<0 ]
Series([], Name: chlHPLC, dtype: float32)
find duplicates
g=chlHPLC2015_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()
doy depth 134 -1.5 3.0 159 -1.5 3.0 Name: depth, dtype: float64
inspect duplicates for doy=134
chlHPLC2015_df[['depth','chlHPLC','doy']]\
[ (chlHPLC2015_df['doy']==134) & (chlHPLC2015_df['depth']==-1.5)]
| depth | chlHPLC | doy | |
|---|---|---|---|
| 10386 | -1.5 | 0.2205 | 134 |
| 10387 | -1.5 | 0.2240 | 134 |
| 10388 | -1.5 | 0.2181 | 134 |
replace with the average
mean=(chlHPLC2015_df[['chlHPLC']] \
[ (chlHPLC2015_df['doy']==134) \
& (chlHPLC2015_df['depth']==-1.5)]).mean()
mean=mean.iloc[0]
chlHPLC2015_df.loc[(chlHPLC2015_df['doy']==134) \
& (chlHPLC2015_df['depth']==-1.5),'chlHPLC']=mean
inspect duplicates for doy=159
chlHPLC2015_df[['depth','chlHPLC','doy']] \
[ (chlHPLC2015_df['doy']==159) & (chlHPLC2015_df['depth']==-1.5)]
| depth | chlHPLC | doy | |
|---|---|---|---|
| 10460 | -1.5 | 0.1700 | 159 |
| 10461 | -1.5 | 0.1684 | 159 |
| 10462 | -1.5 | 0.4612 | 159 |
replace with the average
mean=(chlHPLC2015_df[['chlHPLC']] \
[ (chlHPLC2015_df['doy']==159) \
& (chlHPLC2015_df['depth']==-1.5)]).mean()
mean=mean.iloc[0]
chlHPLC2015_df.loc[(chlHPLC2015_df['doy']==159) \
& (chlHPLC2015_df['depth']==-1.5),'chlHPLC']=mean
chlHPLC2015_df.drop_duplicates(subset=['depth','doy'],inplace=True)
Identify days of year for which the number of observations is less than or equal to 2. There are not measurements at enough depths to be able to calculate a meaningful vertical integration.
g=chlHPLC2015_df.groupby('doy')['doy'].value_counts()
g.where(g<=2).dropna()
doy doy 91 91 2.0 93 93 1.0 96 96 1.0 99 99 1.0 102 102 2.0 106 106 2.0 Name: doy, dtype: float64
chlHPLC2015_df[['depth','chlHPLC','doy']] \
[(chlHPLC2015_df.doy==91) | (chlHPLC2015_df.doy==93) \
| (chlHPLC2015_df.doy==96) | (chlHPLC2015_df.doy==99) \
| (chlHPLC2015_df.doy==102) | (chlHPLC2015_df.doy==106)]
| depth | chlHPLC | doy | |
|---|---|---|---|
| 10306 | -20.0 | 0.0281 | 91 |
| 10307 | -0.5 | 0.0556 | 91 |
| 10308 | -5.0 | 0.0626 | 93 |
| 10309 | -5.0 | 0.0656 | 96 |
| 10310 | -5.0 | 0.0618 | 99 |
| 10311 | -20.0 | 0.0620 | 102 |
| 10312 | -5.0 | 0.1085 | 102 |
| 10316 | -20.0 | 0.0532 | 106 |
| 10317 | -5.0 | 0.0807 | 106 |
By human judgement, I drop these observations.
chlHPLC2015_df=chlHPLC2015_df.drop(
chlHPLC2015_df[(chlHPLC2015_df.doy==91) | (chlHPLC2015_df.doy==93) \
| (chlHPLC2015_df.doy==96) | (chlHPLC2015_df.doy==99) \
| (chlHPLC2015_df.doy==102) | (chlHPLC2015_df.doy==106)].index)
in $mg\ Chla\ m^{-2}$
#doy: 104
depths_case1 =np.array([-5,-20,-65])
weights_case1=np.array([12.5,30,57.5])
#doy: 108
depths_case2 =np.array([-1.5,-5,-20])
weights_case2=np.array([3.25,9.25,87.5])
#doy: 110,112,114,116,159,161
depths_case3 =np.array([-1.5,-5,-20,-40,-60])
weights_case3=np.array([3.25,9.25,17.5,20,50])
#doy: 118
depths_case4 =np.array([-0.5,-1.5,-5,-14,-40,-60])
weights_case4=np.array([1,2.25,9.25,17.5,20,50])
#doy: 120,124,126,128,132,134,136,138,140,142,145,147,149,151,153,155,157
depths_case5 =np.array([-0.5,-1.5,-5,-20,-40,-60])
weights_case5=np.array([1,2.25,9.25,17.5,20,50])
#doy: 122
depths_case6 =np.array([-0.5,-1.5,-5,-40,-60])
weights_case6=np.array([1,2.25,19.25,27.5,50])
#doy: 163
depths_case7 =np.array([-0.5,-1.5,-5,-20,-40])
weights_case7=np.array([1,2.25,9.25,17.5,70])
#doy: 166,169
depths_case8 =np.array([-0.5,-5,-10,-20,-40])
weights_case8=np.array([2.75,4.75,7.5,15,70])
#doy: 172,175,178,181,184,189
depths_case9 =np.array([-0.5,-1.5,-5,-10,-20,-40])
weights_case9=np.array([1,2.25,4.25,7.5,15,70])
#doy: 187
depths_case10 =np.array([-0.5,-1.5,-5,-10,-20,-40,-60])
weights_case10=np.array([1,2.25,4.25,7.5,15,20,50])
#doy: 191
depths_case11 =np.array([-1.5,-5,-8,-18,-40])
weights_case11=np.array([3.25,3.25,6.5,16,71])
#doy: 193
depths_case12 =np.array([-1.5,-15,-30])
weights_case12=np.array([8.25,14.25,77.5])
#doy:195
depths_case13 =np.array([-1.5,-5,-30])
weights_case13=np.array([3.25,14.25,82.5])
nT=365
array1d_iT_obsvintchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
obsvintchlHPLC=np.NaN
chlHPLC=[]
doy=iT+1
chlHPLC_oneday_df=chlHPLC2015_df[['depth','chlHPLC']] \
[(chlHPLC2015_df.doy==doy)]
depths=chlHPLC_oneday_df['depth']
if set(depths_case1)==set(depths):
for depth in depths_case1:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
elif set(depths_case2)==set(depths):
for depth in depths_case2:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
elif set(depths_case3)==set(depths):
for depth in depths_case3:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
elif set(depths_case4)==set(depths):
for depth in depths_case4:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
elif set(depths_case5)==set(depths):
for depth in depths_case5:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
elif set(depths_case6)==set(depths):
for depth in depths_case6:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
elif set(depths_case7)==set(depths):
for depth in depths_case7:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
elif set(depths_case8)==set(depths):
for depth in depths_case8:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
elif set(depths_case9)==set(depths):
for depth in depths_case9:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case9)
elif set(depths_case10)==set(depths):
for depth in depths_case10:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case10)
elif set(depths_case11)==set(depths):
for depth in depths_case11:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case11)
elif set(depths_case12)==set(depths):
for depth in depths_case12:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case12)
elif set(depths_case13)==set(depths):
for depth in depths_case13:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case13)
array1d_iT_obsvintchlHPLC[iT]=obsvintchlHPLC
The bloom initiation was identified when the accumulation rates ($d^{-1}$) were positive for one week (Boss and Behrenfeld, 2010):
accumulation_rate=$\frac{chl_{iT}-chl_{iT-1}}{chl_{iT-1}}$
$In\ situ$ data was smoothed with LOWESS (locally weighted scatterplot smoothing) and interpolated before calculating the accumulation rates.
For $in\ situ$ data, by applying the definition of Boss & Behrenfeld (2010) of 7 days with a positive accumulation rate ($d^{-1}$), the date of bloom initiation would have been April 21. This date was questionable because there was no continuous increase of Chl $a$. This was an indication that human judgment was needed to identify the date of the bloom initiation. We supposed that the small peak observable on April 18 and even the larger one observed on May 12 may not be due to local phytoplankton but instead to the flushing of sea ice algae or to the horizontal advection of phytoplankton. The algorithm was applied after May 12 and found a date of bloom initiation on June 7. Note that Oziel et al. (2019) in $Elem.\ Sci\ Anth.$ also found a date of bloom initiation at the start of June from the same observations.
x=first_year
y=array1d_iT_obsvintchlHPLC
df=pd.DataFrame({'Day of year':x,'Observed Chl a vert. integrated (0-100m)':y})
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,251])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
Figure. Scatterplot of day of year and observed Chl a vert. integrated (0-100m).
y_hat=statsmodels.nonparametric.smoothers_lowess.lowess(y,x,frac=1/4)
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Add the prediction line
fig.add_traces(go.Scatter(x=y_hat[:,0], y=y_hat[:,1], name='LOWESS, frac=1/4', line=dict(color='red')))
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,251])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
# ------- Define interploation functions -------
# Linear - draws a line between the two nearest points and calculates y value based on the slope of that line
f_linear = scipy.interpolate.interp1d(y_hat[:,0],
y=y_hat[:,1],
bounds_error=False)
# Create a new set of points with x values
indices_not_nan=np.argwhere(~np.isnan(y)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
xnew=np.empty(end-start+1)
for iT in range(start,end+1):
xnew[iT-start]=iT
# Find y values based on the interpolation method
array1d_iT_obsvintchlHPLClowess=f_linear(xnew)
start_edge=np.empty(start)
start_edge[:]=np.NaN
end_edge=np.empty(365-1-end)
end_edge[:]=np.NaN
array1d_iT_obsvintchlHPLClowess=np.concatenate((start_edge,
array1d_iT_obsvintchlHPLClowess,
end_edge))
x=first_year
y=array1d_iT_obsvintchlHPLClowess
df=pd.DataFrame({'Day of year':x,'Smoothed and interpolated observed Chl a vert. integrated (0-100m)':y})
indices_not_nan=np.argwhere(~np.isnan(array1d_iT_obsvintchlHPLClowess)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
# no NaN between first and last smoothed observed value
assert np.array_equal(indices_not_nan,np.arange(start,end+1)), (
'There is at least one NaN in array1d_iT_obsvintchlHPLCsmooth before '
'its first and last valid values.')
array1d_iT_r2015=np.empty(365)
array1d_iT_r2015[:]=np.NaN
for iT in range(start+1,end):
chliT=array1d_iT_obsvintchlHPLClowess[iT]
chljT=array1d_iT_obsvintchlHPLClowess[iT-1]
acc=(chliT-chljT)/chljT
array1d_iT_r2015[iT]=acc
x=first_year
y=array1d_iT_r2015
df=pd.DataFrame({'Day of year':x,'Accumulation rates':y})
found=False
iT=138 # starting after the small peak of 131 (May 12) + 7
while not found and iT<=364:
array1d_iT_accweekbefore=array1d_iT_r2015[np.arange(iT-6,iT+1)]
if(np.all(~np.isnan(array1d_iT_accweekbefore)) \
and \
np.all(array1d_iT_accweekbefore>0)):
found=True
else:
iT=iT+1
tI=iT
year=2015 # any non-leap year
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tIs=datetime.datetime(year,1,1) + datetime.timedelta(tI)
Day of year
tI=tI+1
tI
158
tIs.strftime("%B %d")
'June 07'
cI=array1d_iT_obsvintchlHPLClowess[tI-1]
(in $mg\ Chl\ m^{-2}$)
cI
6.280748298393588
blooms_df.loc[2015,['tI']]=tI
blooms_df.loc[2015,['cI']]=cI
with plt.style.context('mplstyles/exp.mplstyle'):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
fig,axs=plt.subplots(2,1,sharex=False,figsize=(16, 6))
# --- ACCUMULATION RATE
ax=axs[0]
h1=ax.plot(first_year,array1d_iT_r2015,'-',color='blue',lw=3)
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks(np.arange(-0.5,1.0,0.5))
ax.set_ylim(-0.5,0.5)
ax.set_ylabel('Accumulation rate ($\mathrm{ d^{-1} }$)')
ax.plot(xlims,[0,0],'k--',color='gray')
# adding the vertical line at start of bloom development
ax.axvline(tI-1,color='green',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[1]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.12,0.15 ,0.65,.60])
axs[1].set_position( [0.12,0.10 ,0.65,.01])
The date of the bloom peak was the date of the maximum (vertically integrated) chlorophyll $a$.
tP=np.nanargmax(array1d_iT_obsvintchlHPLC)
tP=tP.item()
year=2015
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tPs=datetime.datetime(year,1,1) + datetime.timedelta(tP)
(in $mg\ Chl\ m^{-2}$)
cP=np.nanmax(array1d_iT_obsvintchlHPLC)
cP
314.64061871171
Day of year
tP=tP+1
tP
195
tPs.strftime("%B %d")
'July 14'
blooms_df.loc[2015,['tP']]=tP
blooms_df.loc[2015,['cP']]=cP
$m$ is the slope between the bloom initiation (I) and the bloom peak (P) in the log scale. A linear regression of the points between I and P was used.
obsvintchl2015_df=pd.DataFrame(data={'iT':first_year,
'obsvintchl':array1d_iT_obsvintchlHPLC})
obsvintchl2015_df=obsvintchl2015_df[obsvintchl2015_df.obsvintchl.notna()]
y=2015
tI=blooms_df.loc[y,['tI']].item()-1 # doy (1-based) -> doy (0-based)
tP=blooms_df.loc[y,['tP']].item()-1 # doy (1-based) -> doy (0-based)
Points after bloom initiation and before bloom peak
initiation2015_df=obsvintchl2015_df[(obsvintchl2015_df.iT>=tI)
&(obsvintchl2015_df.iT<=tP)]
Points before bloom initiation and after bloom peak. Those points were note used for the linear regression.
out2015_df=obsvintchl2015_df[(obsvintchl2015_df.iT<tI)
|(obsvintchl2015_df.iT>tP)]
loginitiation2015_df=pd.DataFrame(
data={'iT':initiation2015_df.iT,
'obsvintlogchl':np.log10(initiation2015_df.obsvintchl)})
olsres=smf.ols('obsvintlogchl~iT',data=loginitiation2015_df).fit()
m=olsres.params.iT
blooms_df.loc[2015,['m']]=m
To draw the line representing $m$.
olsx2015=[tI,tP]
olsy2015=olsres.predict(pd.DataFrame({'iT':olsx2015}))
olsy2015=np.power(10,olsy2015)
Break-up was on 18 July 2016 (Oziel et al., 2019 in Elem. Sci. Anth.).
obsice=np.empty(365)
obsice[:]=np.NaN
obsice[0:199]=1
obsice[199]=0
Day of in situ sea ice break-up:
obs_breakup=next(i for i,v in enumerate(obsice) if v < 0.50)
The file Nutrients from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ nutrients at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variables 'no3_um_l', 'po4_um_l', and 'sioh4_um_l' are the nitrate, phosphate, and silicic acid concentrations, respectively ($\mathrm{ \mu }$M).
def load_nutrients():
fname='data/DataS2_observations_Qikiqtarjuaq/66412.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'date': str,
'date_time': str,
'latitude': np.float32,
'longitude': np.float32,
'sample_type': 'category',
'sample_source': 'category',
'snow_thickness': 'category',
'code_operation': 'category',
'station': 'category',
'station_type': 'category',
'cast': np.float32,
'bottle': np.float32,
'bottom_depth_m': np.float32,
'target_depth_m': 'string',
'depth_m': np.float32,
'filter_type': 'category',
'quality_flag': 'category',
'dilution_factor': np.float32,
'no3_um_l': np.float32,
'no2_um_l': np.float32,
'no2_and_no3': np.float32,
'po4_um_l': np.float32,
'sioh4_um_l': np.float32,
'n_ice_samples_ice_camp_2016': np.float32,
'method': 'category',
'pi': 'category'
},
parse_dates=['date','date_time']
)
.rename(
columns=dict(
depth_m='depth',
no3_um_l='no3',
po4_um_l='po4',
sioh4_um_l='sioh4'
)
)
)
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
return df
nutrients_df=load_nutrients()
We select the observations of nitrate at the ice camp 2016 in the water column made with a GFF filter.
nutrients_df=nutrients_df[(nutrients_df.mission=='ice_camp_2016') \
& (nutrients_df.sample_type=='water') \
& (nutrients_df.filter_type=='gff')]
Find duplicates.
g=nutrients_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()
Series([], Name: depth, dtype: float64)
No duplicates.
The file Pigments, nutrients, Chlorophyll a and Phaeopigments (concentration) from https://www.seanoe.org/data/00487/59892/ contains $in\ situ$ chlorophyll a at the Qikiqtarjuaq sea ice camps 2015 and 2016 ($67.48^\circ N$, $-63.79^\circ E$). Variables 'conc_mg_m3' is the chlorophyll a concentration measured by HPLC (mg Chl $\mathrm{ m^{-3} }$).
def load_pigment():
fname='data/DataS2_observations_Qikiqtarjuaq/66417.csv'
df = (
pd.read_csv(
fname,
dtype={
'mission': 'string',
'date': str,
'date_time': str,
'latitude': np.float32,
'longitude': np.float32,
'sample_id': 'string',
'sample_code': 'string',
'sample_type': 'category',
'sample_source': 'category',
'snow_thickness': 'category',
'leg': 'category',
'code_operation': 'category',
'station': 'category',
'station_type': 'category',
'cast': np.float32,
'bottle': np.float32,
'target_depth_m': 'string',
'depth_m': np.float32,
'pigment': 'category',
'conc_mg_m2': np.float32,
'conc_ug_l': np.float32,
'conc_mg_m3': np.float32,
'filtered_vol_l': np.float32,
'vol_after_melt_l': np.float32,
'ice_added_vol_of_water': np.float32,
'ice_total_vol_water': np.float32,
'ice_calculated_vol_ice': np.float32,
'dilution_factor': np.float32,
'quality_control_peaks': 'category',
'nb_core': np.float32,
'nb_core_9cm_dia': np.float32,
'core_area_9cm_dia': np.float32,
'nb_core_14cm_dia': np.float32,
'core_area_14cm_dia': np.float32,
'area_m2': np.float32,
'description_of_sample': 'category',
'n': np.float32,
'observations': 'category',
'method': 'category',
'pi': 'category'
},
parse_dates=['date','date_time']
)
.rename(
columns=dict(
depth_m='depth',
conc_ug_l='chlfluo',
conc_mg_m3='chlHPLC'
)
)
)
df['doy']=df['date'].apply(lambda x:x.timetuple().tm_yday)
return df
chlHPLC_df=load_pigment()
We select the observations of chlorophyll a concentrations at the ice camp 2016 measured by HPLC.
chlHPLC_df=chlHPLC_df[(chlHPLC_df.mission=='ice_camp_2016') \
& (chlHPLC_df.sample_type=='water') \
& (chlHPLC_df.pigment=='Total Chlorophyll a') \
& (chlHPLC_df.method=='HPLC')]
remove negative chlorophyll a
chlHPLC_df['chlHPLC'][ chlHPLC_df['chlHPLC']<0 ]
Series([], Name: chlHPLC, dtype: float32)
find duplicates
g=chlHPLC_df.groupby('doy')['depth'].value_counts()
g.where(g>1).dropna()
doy depth
146 1.5 2.0
10.0 2.0
Name: depth, dtype: float64
inspect duplicates for doy=146
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==146) | (chlHPLC_df.doy==148)]
chlHPLC_oneday_df
| depth | chlHPLC | doy | |
|---|---|---|---|
| 20357 | 0.5 | 0.1056 | 146 |
| 20358 | 1.5 | 0.0899 | 146 |
| 20359 | 5.0 | 0.0906 | 146 |
| 20360 | 10.0 | 0.0905 | 146 |
| 20361 | 20.0 | 0.1011 | 146 |
| 20362 | 40.0 | 0.0633 | 146 |
| 20363 | 0.5 | 0.0971 | 148 |
| 20364 | 1.5 | 0.1483 | 146 |
| 20365 | 5.0 | 0.1086 | 148 |
| 20366 | 10.0 | 0.0890 | 146 |
| 20367 | 20.0 | 0.0735 | 148 |
| 20368 | 40.0 | 0.0399 | 148 |
By human judgement, I replace the doy of rows 20364 and 20366 with 148.
chlHPLC_df.loc[[20364,20366],'doy']=148
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==146) | (chlHPLC_df.doy==148)]
chlHPLC_oneday_df
| depth | chlHPLC | doy | |
|---|---|---|---|
| 20357 | 0.5 | 0.1056 | 146 |
| 20358 | 1.5 | 0.0899 | 146 |
| 20359 | 5.0 | 0.0906 | 146 |
| 20360 | 10.0 | 0.0905 | 146 |
| 20361 | 20.0 | 0.1011 | 146 |
| 20362 | 40.0 | 0.0633 | 146 |
| 20363 | 0.5 | 0.0971 | 148 |
| 20364 | 1.5 | 0.1483 | 148 |
| 20365 | 5.0 | 0.1086 | 148 |
| 20366 | 10.0 | 0.0890 | 148 |
| 20367 | 20.0 | 0.0735 | 148 |
| 20368 | 40.0 | 0.0399 | 148 |
Identify days of year for which the number of observations is less than or equal to 2. There are not measurements at enough depths to be able to calculate a meaningful vertical integration.
g=chlHPLC_df.groupby('doy')['doy'].value_counts()
g.where(g<=2).dropna()
doy doy 120 120 2.0 Name: doy, dtype: float64
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC','doy']][(chlHPLC_df.doy==120)]
chlHPLC_oneday_df
| depth | chlHPLC | doy | |
|---|---|---|---|
| 20299 | 1.5 | 0.0355 | 120 |
| 20300 | 5.0 | 0.0253 | 120 |
By human judgement, I drop these observations.
chlHPLC_df=chlHPLC_df.drop(chlHPLC_df[(chlHPLC_df.doy==120)].index)
in $\mathrm{ mg\ Chl\ m^{-2} }$
#doy: 123
depths_case1 =np.array([ 1.5, 5, 10, 20, 40])
weights_case1=np.array([ 3.25,4.25, 7.5, 15, 10])
#doys: 125,127
depths_case2 =np.array([0.5, 1.5, 5, 20, 40])
weights_case2=np.array([1, 2.25,9.25, 17.5, 10])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([0.5,1.5, 5, 10, 20, 40])
weights_case3=np.array([1, 2.25, 4.25, 7.5, 15, 10])
#doy: 141
depths_case4 =np.array([0.5,1.5, 5, 10, 40])
weights_case4=np.array([1, 2.25, 4.25,17.5, 15])
#doy: 174
depths_case5 =np.array([0.5,1.5, 10, 20, 40])
weights_case5=np.array([1, 4.75, 9.25,15, 10])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([0.5,1.5, 10, 20,30, 60])
weights_case6=np.array([1, 4.75, 9.25,10,15, 0])
#doys: 193,195,200
depths_case7 =np.array([ 1.5, 10, 20,30, 60, 75])
weights_case7=np.array([ 5.75, 9.25,10,15, 0, 0])
#doys: 197,204
depths_case8 =np.array([ 1.5, 10, 20,30, 45, 60])
weights_case8=np.array([ 5.75, 9.25,10,12.5, 2.5,0])
nT=365
array1d_iT_obsvint40mchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
obsvintchlHPLC=np.NaN
chlHPLC=[]
doy=iT+1
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
depths=chlHPLC_oneday_df['depth']
if set(depths_case1).issubset(set(depths)):
for depth in depths_case1:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
elif set(depths_case2).issubset(set(depths)):
for depth in depths_case2:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
elif set(depths_case3).issubset(set(depths)):
for depth in depths_case3:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
elif set(depths_case4).issubset(set(depths)):
for depth in depths_case4:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
elif set(depths_case5).issubset(set(depths)):
for depth in depths_case5:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
elif set(depths_case6).issubset(set(depths)):
for depth in depths_case6:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
elif set(depths_case7).issubset(set(depths)):
for depth in depths_case7:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
elif set(depths_case8).issubset(set(depths)):
for depth in depths_case8:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
array1d_iT_obsvint40mchlHPLC[iT]=obsvintchlHPLC
in $\mathrm{ mg\ Chla\ m^{-2} }$
#doy: 123
depths_case1 =np.array([ 1.5, 5, 10, 20, 40])
weights_case1=np.array([ 3.25,4.25, 7.5, 15, 70])
#doys: 125,127
depths_case2 =np.array([0.5, 1.5, 5, 20, 40])
weights_case2=np.array([1, 2.25,9.25, 17.5, 70])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([0.5,1.5, 5, 10, 20, 40])
weights_case3=np.array([1, 2.25, 4.25, 7.5, 15, 70])
#doy: 141
depths_case4 =np.array([0.5,1.5, 5, 10, 40])
weights_case4=np.array([1, 2.25, 4.25,17.5, 75])
#doy: 174
depths_case5 =np.array([0.5,1.5, 10, 20, 40])
weights_case5=np.array([1, 4.75, 9.25,15, 70])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([0.5,1.5, 10, 20,30, 60])
weights_case6=np.array([1, 4.75, 9.25,10,20, 55])
#doys: 193,195,200
depths_case7 =np.array([ 1.5, 10, 20,30, 60, 75])
weights_case7=np.array([ 5.75, 9.25,10,20, 22.5,32.5])
#doys: 197,204
depths_case8 =np.array([ 1.5, 10, 20,30, 45,60])
weights_case8=np.array([ 5.75, 9.25,10,12.5,15,47.5])
nT=365
array1d_iT_obsvintchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
obsvintchlHPLC=np.NaN
chlHPLC=[]
doy=iT+1
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
depths=chlHPLC_oneday_df['depth']
if set(depths_case1).issubset(set(depths)):
for depth in depths_case1:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case1)
elif set(depths_case2).issubset(set(depths)):
for depth in depths_case2:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case2)
elif set(depths_case3).issubset(set(depths)):
for depth in depths_case3:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case3)
elif set(depths_case4).issubset(set(depths)):
for depth in depths_case4:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case4)
elif set(depths_case5).issubset(set(depths)):
for depth in depths_case5:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case5)
elif set(depths_case6).issubset(set(depths)):
for depth in depths_case6:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case6)
elif set(depths_case7).issubset(set(depths)):
for depth in depths_case7:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case7)
elif set(depths_case8).issubset(set(depths)):
for depth in depths_case8:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsvintchlHPLC=np.dot(chlHPLC,weights_case8)
array1d_iT_obsvintchlHPLC[iT]=obsvintchlHPLC
in $mg\ Chla\ m^{-3}$
#doy: 123
depths_case1 =np.array([ 1.5, 5, 10, 20, 40])
weights_case1=np.array([ 3.25, 4.25, 7.5, 15, 10])
#doys: 125,127
depths_case2 =np.array([ 0.5, 1.5, 5, 20, 40])
weights_case2=np.array([ 1, 2.25, 9.25, 17.5, 10])
#doys:130,132,134,137,139,144,146,148,151,153,155,158,160,162,165,167,169,172,176
depths_case3 =np.array([ 0.5, 1.5, 5, 10, 20, 40])
weights_case3=np.array([ 1, 2.25, 4.25, 7.5, 15, 10])
#doy: 141
depths_case4 =np.array([ 0.5, 1.5, 5, 10, 40])
weights_case4=np.array([ 1, 2.25, 4.25, 17.5, 15])
#doy: 174
depths_case5 =np.array([ 0.5, 1.5, 10, 20, 40])
weights_case5=np.array([ 1, 4.75, 9.25, 15, 10])
#doys: 179,181,183,186,188,190
depths_case6 =np.array([ 0.5, 1.5, 10, 20, 30, 60])
weights_case6=np.array([ 1, 4.75, 9.25, 10, 15, 0])
#doys: 193,195,200
depths_case7 =np.array([ 1.5, 10, 20, 30, 60, 75])
weights_case7=np.array([ 5.75, 9.25, 10, 15, 0, 0])
#doys: 197,204
depths_case8 =np.array([ 1.5, 10, 20, 30, 45, 60])
weights_case8=np.array([ 5.75, 9.25, 10, 12.5, 2.5, 0])
nT=365
array1d_iT_obsmeanchlHPLC=np.full(nT,np.NaN)
for iT in np.arange(0,nT):
obsmeanchlHPLC=np.NaN
chlHPLC=[]
doy=iT+1
chlHPLC_oneday_df=chlHPLC_df[['depth','chlHPLC']][(chlHPLC_df.doy==doy)]
depths=chlHPLC_oneday_df['depth']
if set(depths_case1).issubset(set(depths)):
for depth in depths_case1:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case1)
elif set(depths_case2).issubset(set(depths)):
for depth in depths_case2:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case2)
elif set(depths_case3).issubset(set(depths)):
for depth in depths_case3:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case3)
elif set(depths_case4).issubset(set(depths)):
for depth in depths_case4:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case4)
elif set(depths_case5).issubset(set(depths)):
for depth in depths_case5:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case5)
elif set(depths_case6).issubset(set(depths)):
for depth in depths_case6:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case6)
elif set(depths_case7).issubset(set(depths)):
for depth in depths_case7:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case7)
elif set(depths_case8).issubset(set(depths)):
for depth in depths_case8:
chlHPLCtempo \
=chlHPLC_oneday_df['chlHPLC'][(chlHPLC_oneday_df.depth==depth)]
assert chlHPLCtempo.size==1
chlHPLCtempo=chlHPLCtempo.iat[0]
chlHPLC.append(chlHPLCtempo)
obsmeanchlHPLC=np.average(np.array(chlHPLC),weights=weights_case8)
array1d_iT_obsmeanchlHPLC[iT]=obsmeanchlHPLC
The bloom initiation was identified when the accumulation rates ($d^{-1}$) were positive for one week (Boss and Behrenfeld, 2010):
accumulation_rate=$\frac{chl_{iT}-chl_{iT-1}}{chl_{iT-1}}$
$In\ situ$ data was smoothed with LOWESS (locally weighted scatterplot smoothing) and interpolated before calculating the accumulation rates.
For $in\ situ$ data, by applying the definition of Boss & Behrenfeld (2010) of 7 days with a positive accumulation rate ($d^{-1}$), the date of bloom initiation would have been May 9. This date was questionable because light was low on May 9, there was no continuous decrease of nutrients and there was no continuous increase of Chl $a$. Also, changing the number of days with a positive accumulation rate from 7 days to 8 days changed the date of bloom initiation from May 9 to May 30 indicating a high sensitivity on the number of days with a positive accumulation rate. This was an indication that human judgment was needed to identify the date of the bloom initiation. We supposed that the small peak observable on May 9 may not be due to local phytoplankton but instead to the flushing of sea ice algae or to the horizontal advection of phytoplankton. The algorithm was applied after May 9 and found a date of bloom initiation on May 30. Note that Oziel et al. (2019) in $Elem.\ Sci\ Anth.$ also found a date of bloom initiation at the end of May from the same observations.
x=first_year
y=array1d_iT_obsvintchlHPLC
df=pd.DataFrame({'Day of year':x,'Observed Chl a vert. integrated (0-100m)':y})
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,251])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
Figure. Scatterplot of day of year and observed Chl a vert. integrated (0-100m).
y_hat=statsmodels.nonparametric.smoothers_lowess.lowess(y,x,frac=1/4)
# Create a scatter plot
fig = px.scatter(df, x=df['Day of year'], y=df['Observed Chl a vert. integrated (0-100m)'],
opacity=0.8, color_discrete_sequence=['black'])
# Add the prediction line
fig.add_traces(go.Scatter(x=y_hat[:,0], y=y_hat[:,1], name='LOWESS, frac=1/4', line=dict(color='red')))
# Change chart background color
fig.update_layout(dict(plot_bgcolor = 'white'))
# Update axes lines
fig.update_xaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black')
fig.update_yaxes(showgrid=True, gridwidth=1, gridcolor='lightgrey',
zeroline=True, zerolinewidth=1, zerolinecolor='lightgrey',
showline=True, linewidth=1, linecolor='black',
range=[0,251])
# Set figure title
fig.update_layout(title=dict(text="Observed Chl a vert. integrated (0-100m)",
font=dict(color='black')))
# Update marker size
fig.update_traces(marker=dict(size=3))
fig.show()
# ------- Define interploation functions -------
# Linear - draws a line between the two nearest points and calculates y value based on the slope of that line
f_linear = scipy.interpolate.interp1d(y_hat[:,0],
y=y_hat[:,1],
bounds_error=False)
# Create a new set of points with x values
indices_not_nan=np.argwhere(~np.isnan(y)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
xnew=np.empty(end-start+1)
for iT in range(start,end+1):
xnew[iT-start]=iT
# Find y values based on the interpolation method
array1d_iT_obsvintchlHPLClowess=f_linear(xnew)
start_edge=np.empty(start)
start_edge[:]=np.NaN
end_edge=np.empty(365-1-end)
end_edge[:]=np.NaN
array1d_iT_obsvintchlHPLClowess=np.concatenate((start_edge,
array1d_iT_obsvintchlHPLClowess,
end_edge))
x=first_year
y=array1d_iT_obsvintchlHPLClowess
df=pd.DataFrame({'Day of year':x,'Smoothed and interpolated observed Chl a vert. integrated (0-100m)':y})
indices_not_nan=np.argwhere(~np.isnan(array1d_iT_obsvintchlHPLClowess)).flatten()
start=indices_not_nan[0]
end =indices_not_nan[-1]
# no NaN between first and last smoothed observed value
assert np.array_equal(indices_not_nan,np.arange(start,end+1)), (
'There is at least one NaN in array1d_iT_obsvintchlHPLCsmooth before '
'its first and last valid values.')
array1d_iT_r2016=np.empty(365)
array1d_iT_r2016[:]=np.NaN
for iT in range(start+1,end):
chliT=array1d_iT_obsvintchlHPLClowess[iT]
chljT=array1d_iT_obsvintchlHPLClowess[iT-1]
acc=(chliT-chljT)/chljT
array1d_iT_r2016[iT]=acc
x=first_year
y=array1d_iT_r2016
df=pd.DataFrame({'Day of year':x,'Accumulation rates':y})
found=False
iT=136 # starting after the small peak of 129 (May 9) + 7
while not found and iT<=364:
array1d_iT_accweekbefore=array1d_iT_r2016[np.arange(iT-6,iT+1)]
if(np.all(~np.isnan(array1d_iT_accweekbefore)) \
and \
np.all(array1d_iT_accweekbefore>0)):
found=True
else:
iT=iT+1
tI=iT
year=2016 # any leap year
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tIs=datetime.datetime(year,1,1) + datetime.timedelta(tI)
Day of year
tI=tI+1
tI
148
tIs.strftime("%B %d")
'May 27'
cI=array1d_iT_obsvintchlHPLClowess[tI-1]
(in $mg\ Chl\ m^{-2}$)
cI
6.412508830445615
blooms_df.loc[2016,['tI']]=tI
blooms_df.loc[2016,['cI']]=cI
with plt.style.context('mplstyles/exp.mplstyle'):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
fig,axs=plt.subplots(2,1,sharex=False,figsize=(16, 6))
# --- ACCUMULATION RATE
ax=axs[0]
h1=ax.plot(first_year,array1d_iT_r2016,'-',color='blue',lw=3)
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks(np.arange(-0.5,1.0,0.5))
ax.set_ylim(-0.5,0.5)
ax.set_ylabel('Accumulation rate ($\mathrm{ d^{-1} }$)')
ax.plot(xlims,[0,0],'k--',color='gray')
# adding the vertical line at start of bloom development
ax.axvline(tI-1,color='green',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[1]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.12,0.15 ,0.65,.60])
axs[1].set_position( [0.12,0.10 ,0.65,.01])
The date of the bloom peak was the date of the maximum (vertically integrated) chlorophyll $a$.
tP=np.nanargmax(array1d_iT_obsvintchlHPLC)
tP=tP.item()
year=2016
# transform day of year into date
# https://stackoverflow.com/questions/2427555/python-question-year-and-day-of-year-to-date
# see answer of Ignacio Vazquez-Abrams
tPs=datetime.datetime(year,1,1) + datetime.timedelta(tP)
(in $mg\ Chl\ m^{-2}$)
cP=np.nanmax(array1d_iT_obsvintchlHPLC)
cP
204.44257855415344
Day of year
tP=tP+1
tP
197
tPs.strftime("%B %d")
'July 15'
blooms_df.loc[2016,['tP']]=tP
blooms_df.loc[2016,['cP']]=cP
$m$ is the slope between the bloom initiation (I) and the bloom peak (P) in the log scale. A linear regression of the points between I and P was used.
obsvintchl2016_df=pd.DataFrame(data={'iT':first_year,
'obsvintchl':array1d_iT_obsvintchlHPLC})
obsvintchl2016_df=obsvintchl2016_df[obsvintchl2016_df.obsvintchl.notna()]
y=2016
tI=blooms_df.loc[y,['tI']].item()-1 # doy (1-based) -> doy (0-based)
tP=blooms_df.loc[y,['tP']].item()-1 # doy (1-based) -> doy (0-based)
Points after bloom initiation and before bloom peak
initiation2016_df=obsvintchl2016_df[(obsvintchl2016_df.iT>=tI)
&(obsvintchl2016_df.iT<=tP)]
Points before bloom initiation and after bloom peak. Those points were note used for the linear regression.
out2016_df=obsvintchl2016_df[(obsvintchl2016_df.iT<tI)
|(obsvintchl2016_df.iT>tP)]
loginitiation2016_df=pd.DataFrame(
data={'iT':initiation2016_df.iT,
'obsvintlogchl':np.log10(initiation2016_df.obsvintchl)})
olsres=smf.ols('obsvintlogchl~iT',data=loginitiation2016_df).fit()
m=olsres.params.iT
blooms_df.loc[2016,['m']]=m
To draw the line representing $m$.
olsx2016=[tI,tP]
olsy2016=olsres.predict(pd.DataFrame({'iT':olsx2016}))
olsy2016=np.power(10,olsy2016)
input file: 1D_GDPS_qsw_y2016.nc
name: qsw
Variable qsw is the shortwave radiation just above surface. It is at indices (2,2) meaning center of grid point was used. File 1D_GDPS_qsw_y2016.nc is the output of CMC GDPS reforecasts (CGRF, Smith et al., 2014, https://doi.org/10.1002/qj.2194.) CMC: Canadian Meteorological Centre's. GDPS: global deterministic prediction system. The data was provided by Gregory Smith (Environment Canada). Gaetan Olivier (UBO) preprocessed it to get the data at the Green Edge 2016 ice camp (67.48N, -63.79E).
The goal was to calculate PAR just below open water ($I_w(0^-)$, $\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$) from shortwave radiation just above surface ($Q(0^+)$, $\mathrm{ W\ m^{-2} }$).
ndays=365
first_year=np.arange(0,ndays)
qswow_file='data/DataS6_output_cgrf/1D_GDPS_qsw_y2016.nc'
array1d_iT_qswowabovehourly=netcdf_tools.read_netcdf(qswow_file,'qsw')
array1d_iT_qswowabovehourly=array1d_iT_qswowabovehourly[:,1,1]
array1d_idoy_qswowabove=np.zeros(ndays)
for idoy in first_year:
qswowabovedaily=0
for iT in np.arange(idoy*24,(idoy+1)*24):
qswowabovedaily=qswowabovedaily+array1d_iT_qswowabovehourly[iT]
qswowabovedaily=qswowabovedaily/24
array1d_idoy_qswowabove[idoy]=qswowabovedaily
albedo=0.066 # model PISCES, ref.: Aumont et al. (2015)
array1d_idoy_qswow=array1d_idoy_qswowabove*(1-albedo)
Wm2_to_PAR=0.43 # model PISCES, ref.: Aumont et al. (2015)
array1d_idoy_PARow=array1d_idoy_qswow*Wm2_to_PAR
onecolor=1/3
array1d_idoy_blueWm2ow=array1d_idoy_PARow*onecolor
array1d_idoy_greenWm2ow=array1d_idoy_PARow*onecolor
array1d_idoy_redWm2ow=array1d_idoy_PARow*onecolor
ref: https://gud.mit.edu/MITgcm/source/pkg/gud/gud_radtrans_direct.F?v=gud
rmus=1/0.83
planck = 6.6256e-34 # Plancks constant J sec
c = 2.998e8 # speed of light m/sec
hc = 1.0/(planck*c)
oavo = 1.0/6.023e23 # 1/Avogadros number
hcoavo = hc*oavo
wb_center=np.array([450,550,650]) # [blue, green, red]
rlamm=wb_center*1e-9 # lambda in m
WtouEins=1e6*rlamm*hcoavo
array1d_idoy_blueow =array1d_idoy_blueWm2ow *WtouEins[0]
array1d_idoy_greenow=array1d_idoy_greenWm2ow*WtouEins[1]
array1d_idoy_redow =array1d_idoy_redWm2ow *WtouEins[2]
array1d_idoy_PARowuEin=array1d_idoy_blueow+array1d_idoy_greenow+array1d_idoy_redow
The file siarea.nemo.2016.365.32bits.bin contains data from a simulation of NEMO3.6 LIM3.6 configured by Gaetan Olivier (UBO) for the Qikiqtarjuaq sea ice camp location ($67.48^\circ N$, $-63.79^\circ E$) in 2016. It contains sea ice concentration (unitless).
ice=np.fromfile('data/DataS4_output_nemo_lim3/siarea.nemo.2016.365.32bits.bin',
dtype='>f')
input file: GE_mod_var.nc.
long_name: solar_heat_flux_under_ice_for_100:100_ice_cover.
standard_name in NEMO: surface_downwelling_shortwave_in water.
Variable solar_heat_flux_under_ice_for_100:100_ice_cover is the surface downwelling shortwave in water (below ice for 100% ice cover). It is at indices (2,2) meaning center of grid point was used. File GE_mod_var.nc is the output of NEMO3.6+LIM3.6 model configured by Gaetan Olivier (UBO) for the Green Edge ice camp 2016 (67.48N, -63.79E).
NEMO3.6: Madec et al., 2017: http://hdl.handle.net/2122/13309.
LIM3.6: Rousset et al., 2015: doi:10.5194/gmdd-8-3403-2015.
The goal was to calculate PAR just below sea ice ($I_i(0^-)$, $\mathrm{ \mu mol\ photons\ m^{-2}\ s^{-1} }$) from shortwave radiation just below sea ice ($Q_i(0^-)$, $\mathrm{ W\ m^{-2} }$).
from 2013 to 2017
qswice_file='data/DataS4_output_nemo_lim3/GE_mod_var.nc'
array1d_iT_qswicefull=netcdf_tools.read_netcdf(qswice_file,'solar heat flux under ice for 100:100 ice cover')
array1d_iT_qswicefull=array1d_iT_qswicefull[:,1,1]
in 2016
array1d_idoy_qswice=array1d_iT_qswicefull[365*3:365*4]
Wm2_to_PAR=0.85 # Lebrun et al. (submitted)
array1d_idoy_PARice=array1d_idoy_qswice*Wm2_to_PAR
# Lebrun et al. (submitted). Values for the Green Edge ice camp 2016.
blue_fraction =0.43
green_fraction=0.46
red_fraction =0.11
array1d_idoy_blueWm2 =array1d_idoy_PARice*blue_fraction
array1d_idoy_greenWm2=array1d_idoy_PARice*green_fraction
array1d_idoy_redWm2 =array1d_idoy_PARice*red_fraction
array1d_idoy_blue =array1d_idoy_blueWm2 *WtouEins[0]
array1d_idoy_green=array1d_idoy_greenWm2*WtouEins[1]
array1d_idoy_red =array1d_idoy_redWm2 *WtouEins[2]
array1d_idoy_PARiceuEin=array1d_idoy_blue+array1d_idoy_green \
+array1d_idoy_red
array1d_idoy_PARiceuEin=array1d_idoy_PARiceuEin.squeeze()
The file siarea.nemo.2016.365.32bits.bin contains data from a simulation of NEMO3.6+LIM3.6 configured by Gaetan Olivier (UBO) for the Qikiqtarjuaq sea ice camp location (67.48N, -63.79E) in 2016. It contains sea ice concentration (unitless).
ice=np.fromfile('data/DataS4_output_nemo_lim3/siarea.nemo.2016.365.32bits.bin',
dtype='>f')
indir='data/DataS8_output_mitgcm/exp0'
info=gudinfo.RunInfo(indir)
$K_{NO_3}=a_{K_{NO_3}}V^{b_{K_{NO_3}}}$
nplk=info.nplk
K=np.empty(nplk)
K[:]=np.nan
a_kn_no3=info.params['a_kn_no3'] # in mmol N m^3
vol=info.f['vol'] # in um^3
b_kn_no3=info.params['b_kn_no3'] # unitless
for iplk in range(0,nplk):
grp=info.f['grp'][iplk]
if (0<=grp and grp<8):
K[iplk]=a_kn_no3[grp]*vol[iplk]**b_kn_no3[grp]
$Q^{min}_{N}=a_{Q^{min}_{N}}V^{b_{Q^{min}_{N}}}$
nplk=info.nplk
Q_min=np.empty(nplk)
Q_min[:]=np.nan
a_qmin_n=info.params['a_qmin_n'] # in mmol N (mmol C)^-1
vol=info.f['vol'] # in um^3
b_qmin_n=info.params['b_qmin_n'] # unitless
for iplk in range(0,nplk):
grp=info.f['grp'][iplk]
if (0<=grp and grp<8):
Q_min[iplk]=a_qmin_n[grp]*vol[iplk]**b_qmin_n[grp]
$V^{max}_{NO_3}=a_{V^{max}_{NO_3}}V^{b_{V^{max}_{NO_3}}}$
nplk=info.nplk
V_max_d=np.empty(nplk)
V_max_d[:]=np.nan
a_vmax_no3_s=info.params['a_vmax_no3'] # in mmol N (mmol C)^-1 s^-1
a_vmax_no3_d=a_vmax_no3_s*86400 # in mmol N (mmol C)^-1 d^-1
vol=info.f['vol'] # in um^3
b_vmax_no3=info.params['b_vmax_no3'] # unitless
for iplk in range(0,nplk):
grp=info.f['grp'][iplk]
if (0<=grp and grp<8):
V_max_d[iplk]=a_vmax_no3_d[grp]*vol[iplk]**b_vmax_no3[grp]
carfile='data/DataS8_output_mitgcm/exp0/car.0000000000.t001.nc'
array2d_idepth_iT_modno3full=netcdf_tools.read_netcdf(carfile, 'TRAC04').squeeze().transpose()
array2d_idepth_iT_modno3=array2d_idepth_iT_modno3full[:,array1d_iT1y_iT]
array2d_idepth_iT_modsioh4full=netcdf_tools.read_netcdf(carfile, 'TRAC06').squeeze().transpose()
array2d_idepth_iT_modsioh4=array2d_idepth_iT_modsioh4full[:,array1d_iT1y_iT]
array2d_idepth_iT_modpo4full=netcdf_tools.read_netcdf(carfile, 'TRAC05').squeeze().transpose()
array2d_idepth_iT_modpo4=array2d_idepth_iT_modpo4full[:,array1d_iT1y_iT]
The values at -0.51 m on the model grid are extrapolated to 0 m. It is the only extrapolation.
array1d_iobs_no3mod=list()
for iobs,row in nutrients_df.iterrows():
doy=row['doy']
depthobs=row['depth']
assert depthobs>=0, \
"depthobs is negative: %r" %(depthobs)
assert depthobs<RC[-1], \
"depthobs %r is below maximal depth of model grid %r" %(depthobs,RC[-1])
y_interp=scipy.interpolate.interp1d(
RC,
array2d_idepth_iT_modno3[:,doy-1].flatten(),
bounds_error=False,
fill_value="extrapolate"
)
no3mod=y_interp(depthobs).item()
array1d_iobs_no3mod.append(no3mod)
no3mod_df=pd.DataFrame({
'depth':nutrients_df['depth'],
'no3mod':array1d_iobs_no3mod,
'doy':nutrients_df['doy']
})
array1d_iobs_sioh4mod=list()
for iobs,row in nutrients_df.iterrows():
doy=row['doy']
depthobs=row['depth']
assert depthobs>=0, \
"depthobs is positive: %r" %(depthobs)
assert depthobs<RC[-1], \
"depthobs %r is below maximal depth of model grid %r" %(depthobs,RC[-1])
y_interp=scipy.interpolate.interp1d(
RC,
array2d_idepth_iT_modsioh4[:,doy-1].flatten(),
bounds_error=False,
fill_value="extrapolate"
)
sioh4mod=y_interp(depthobs).item()
array1d_iobs_sioh4mod.append(sioh4mod)
sioh4mod_df=pd.DataFrame({
'depth':nutrients_df['depth'],
'sioh4mod':array1d_iobs_sioh4mod,
'doy':nutrients_df['doy']
})
array1d_iobs_po4mod=list()
for iobs,row in nutrients_df.iterrows():
doy=row['doy']
depthobs=row['depth']
assert depthobs>=0, \
"depthobs is positive: %r" %(depthobs)
assert depthobs<RC[-1], \
"depthobs %r is below maximal depth of model grid %r" %(depthobs,RC[-1])
y_interp=scipy.interpolate.interp1d(
RC,
array2d_idepth_iT_modpo4[:,doy-1].flatten(),
bounds_error=False,
fill_value="extrapolate"
)
po4mod=y_interp(depthobs).item()
array1d_iobs_po4mod.append(po4mod)
po4mod_df=pd.DataFrame({
'depth':nutrients_df['depth'],
'po4mod':array1d_iobs_po4mod,
'doy':nutrients_df['doy']
})
chlfile='data/DataS8_output_mitgcm/exp0/chl.0000000000.t001.nc'
array2d_idepth_iT_modchlfull=\
read_mitgcm.get_array2d_idepth_iT_chlfull(chlfile)
array2d_idepth_iT_modchl=array2d_idepth_iT_modchlfull[:,array1d_iT1y_iT]
array1d_iobs_chlmod=list()
for iobs,row in chlHPLC_df.iterrows():
doy=row['doy']
depthobs=row['depth']
assert depthobs>0, \
"depthobs is not negative: %r" %(depthobs)
assert depthobs<RC[-1], \
"depthobs %r is below maximal depth of model grid %r" %(depthobs,RC[-1])
y_interp=scipy.interpolate.interp1d(
RC,
array2d_idepth_iT_modchl[:,doy-1].flatten(),
bounds_error=False,
fill_value="extrapolate"
)
chlmod=y_interp(depthobs).item()
array1d_iobs_chlmod.append(chlmod)
chlmod_df=pd.DataFrame({
'depth':chlHPLC_df['depth'],
'chlmod':array1d_iobs_chlmod,
'doy':chlHPLC_df['doy']
})
in $\mathrm{ mg\ Chla\ m^{-2} }$
array1d_iT_modvint40mchl=vstats.vint(
array2d_idepth_iT_tracer=array2d_idepth_iT_modchl,
array1d_idepth_delR=drF,
depth_end=40)
in $\mathrm{ mg\ Chla\ m^{-2} }$
array1d_iT_modvintchl=vstats.vint(
array2d_idepth_iT_tracer=array2d_idepth_iT_modchl,
array1d_idepth_delR=drF,
depth_end=100)
in $\mathrm{ mg\ Chla\ m^{-3} }$
array1d_iT_modmeanchl=vstats.mean(
array2d_idepth_iT_tracer=array2d_idepth_iT_modchl,
array1d_idepth_delR=drF,
depth_end=40)
biofile='data/DataS8_output_mitgcm/exp0/car.0000000000.t001.nc'
array2d_idepth_iT_modphytobiofull=\
read_mitgcm.get_array2d_idepth_iT_phytobiofull(biofile)
array2d_idepth_iT_modphytobio=\
array2d_idepth_iT_modphytobiofull[:,array1d_iT1y_iT]
in $\mathrm{ mg\ C\ m^{-2} }$
array1d_iT_modvint40mphytobio=vstats.vint(
array2d_idepth_iT_tracer=array2d_idepth_iT_modphytobio,
array1d_idepth_delR=drF,
depth_end=40)
in $\mathrm{ mg\ C\ m^{-2} }$
array1d_iT_modvintphytobio=vstats.vint(
array2d_idepth_iT_tracer=array2d_idepth_iT_modphytobio,
array1d_idepth_delR=drF,
depth_end=100)
in $\mathrm{ mg\ C\ m^{-3} }$
array1d_iT_modmeanphytobio=vstats.mean(
array2d_idepth_iT_tracer=array2d_idepth_iT_modphytobio,
array1d_idepth_delR=drF,
depth_end=40)
type_id=np.arange(1,52)
pos_Prochlorococcus=np.arange( 0, 1)
pos_Synechococcus =np.arange( 1, 2)
pos_smalleuk =np.arange( 2, 4)
pos_cocco =np.arange( 4, 9)
pos_diazo =np.arange( 9,13)
pos_Trichodesmium =np.arange(13,14)
pos_diatoms =np.arange(14,23)
pos_mixodino =np.arange(23,33)
pos_zoo =np.arange(33,49)
def make_plots(axs):
box=dict(pad=5,alpha=0)
# --- CELL NUTRIENT UPTAKE HALF SATURATION CONSTANT ON NITRATE
ax=axs[0]
h3=ax.plot(type_id[pos_smalleuk]-2, K[pos_smalleuk],'gh',
markersize=12)
h4=ax.plot(type_id[pos_cocco]-2, K[pos_cocco], 'b+')
h7=ax.plot(type_id[pos_diatoms]-2-5, K[pos_diatoms], 'ro')
h8=ax.plot(type_id[pos_mixodino]-2-5,K[pos_mixodino],'md')
h9=ax.plot(type_id[pos_zoo]-2-5, K[pos_zoo], 'k*')
ax.set_xticklabels([])
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator( 1))
ax.set_xlim([0,43])
ax.set_ylabel('$K_{NO_3}\ (\mathrm{ mmol\ N\ m^{-3} })$')
ax.set_yscale('log')
ax.set_ylim([0.1,20])
ax.grid()
plt.text(0.05,0.92,'(a)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- CELL MINIMUM STOICHIOMETRIC QUOTAS OF NITROGEN RELATIVE TO CARBON
ax=axs[1]
h3=ax.plot(type_id[pos_smalleuk]-2, Q_min[pos_smalleuk],'gh',
markersize=12)
h4=ax.plot(type_id[pos_cocco]-2, Q_min[pos_cocco], 'b+')
h7=ax.plot(type_id[pos_diatoms]-2-5, Q_min[pos_diatoms], 'ro')
h8=ax.plot(type_id[pos_mixodino]-2-5,Q_min[pos_mixodino],'md')
h9=ax.plot(type_id[pos_zoo]-2-5, Q_min[pos_zoo], 'k*')
ax.set_xticklabels([])
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator( 1))
ax.set_xlim([0,43])
ax.set_ylabel('$Q^{min}_N\ (\mathrm{ mmol\ N\ (mmol\ C)^{-1} })$')
ax.set_yscale('log')
ax.set_ylim([0.004,0.1])
ax.grid()
plt.text(0.05,0.92,'(b)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- CELL NUTRIENT UPTAKE RATE OF NITRATE RELATIVE TO CARBON
ax=axs[2]
h3=ax.plot(type_id[pos_smalleuk]-2, V_max_d[pos_smalleuk],'gh',
markersize=12,label='picophytoplankton')
h4=ax.plot(type_id[pos_cocco]-2, V_max_d[pos_cocco], 'b+',
label='other nanophytoplankton')
h7=ax.plot(type_id[pos_diatoms]-2-5, V_max_d[pos_diatoms], 'ro',
label='diatoms')
h8=ax.plot(type_id[pos_mixodino]-2-5,V_max_d[pos_mixodino],'md',
label='mixotrophic dinoflagellates')
h9=ax.plot(type_id[pos_zoo]-2-5, V_max_d[pos_zoo], 'k*',
label='microzooplankton')
ax.set_xlabel('Planktonic type ID',fontsize=23)
ax.xaxis.set_major_locator(mticker.MultipleLocator(10))
ax.xaxis.set_minor_locator(mticker.MultipleLocator( 1))
ax.set_xlim([0,43])
ax.set_ylabel('$V^{max}_{NO_3}\ (\mathrm{ mmol\ N\ (mmol\ C)^{-1}\ d^{-1} })$')
ax.set_yscale('log')
ax.set_ylim([0.005,1])
ax.grid()
ax.legend(loc='lower right')
plt.text(0.05,0.92,'(c)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# --- POSITION
axs[0].set_position( [0.10,0.69 ,0.85,0.29])
axs[1].set_position( [0.10,0.38 ,0.85,0.29])
axs[2].set_position( [0.10,0.07 ,0.85,0.29])
with plt.style.context('mplstyles/gud_traits.mplstyle'):
# Plot
fig,axs=plt.subplots(3,1,figsize=(16, 11.6))
make_plots(axs)
fig.align_ylabels(axs[:])
# --- SAVE
plt.savefig('figures_supmat_progress/traits_supp.png',dpi=1000)
def make_plots(axs):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
# --- DAILY SHORTWAVE RADIATION JUST ABOVE SURFACE IN 2016
ax=axs[0,0]
ax.plot(first_year,array1d_idoy_qswowabove,color='black',lw=3)
ax.set_ylabel('$Q$($\mathrm{0^+}$)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,400)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.05,
'(a) shortwave radiation just above surface',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- BLANK
ax=axs[0,1]
ax.axis('off')
# --- DAILY SHORTWAVE RADIATION JUST BELOW SEA ICE FROM 2016
ax=axs[1,0]
ax.plot(first_year,array1d_idoy_qswice,color='black',lw=3)
ax.set_ylabel('$Q_\mathrm{i}$($\mathrm{0^-}$)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,40)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.40,
'JUST BELOW SEA ICE',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.50,1.10,
'(b) shortwave radiation',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY SHORTWAVE RADIATION JUST BELOW OPEN WATER IN 2016
ax=axs[1,1]
ax.plot(first_year,array1d_idoy_qswow,color='black',lw=3)
ax.set_ylabel('$Q_\mathrm{w}$($\mathrm{0^-}$)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,400)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.40,
'JUST BELOW OPEN WATER',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
plt.text(0.50,1.10,
'(c) shortwave radiation',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY PAR JUST BELOW SEA ICE IN 2016
# --- W m^-2
ax=axs[2,0]
ax.plot(first_year,array1d_idoy_PARice,color='black',lw=3)
ax.set_ylabel('$I_\mathrm{i}$($\mathrm{0^-}$)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,40)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(d) PAR',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY PAR JUST BELOW OPEN WATER IN 2016
# --- W m^-2
ax=axs[2,1]
ax.plot(first_year,array1d_idoy_PARow,color='black',lw=3)
ax.set_ylabel('$I_\mathrm{w}$($\mathrm{0^-}$)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,400)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(e) PAR',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY IRRADIANCE OF EACH 3 RGB BANDS JUST BELOW SEA ICE IN 2016
# --- W m^-2
ax=axs[3,0]
ax.plot(first_year,array1d_idoy_blueWm2,
':',color='blue',lw=3,label='blue')
ax.plot(first_year,array1d_idoy_greenWm2,
'--',color='green',lw=3,label='green')
ax.plot(first_year,array1d_idoy_redWm2,
'-.',color='red',lw=3,label='red')
ax.legend()
ax.set_ylabel('$E_\mathrm{i}$($\mathrm{0^-}$)(band)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,40)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(f) irradiance of each RGB band',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY IRRADIANCE OF EACH 3 RGB BANDS JUST BELOW OPEN WATER IN 2016
# --- W m^-2
ax=axs[3,1]
ax.plot(first_year,array1d_idoy_blueWm2ow,
':',color='blue',lw=3,label='blue')
ax.plot(first_year,array1d_idoy_greenWm2ow,
'--',color='green',lw=3,label='green')
ax.plot(first_year,array1d_idoy_redWm2ow,
'-.',color='red',lw=3,label='red')
ax.legend()
ax.set_ylabel('$E_\mathrm{w}$($\mathrm{0^-}$)(band)\n'
'(W $\mathrm{m^{-2}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,400)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(g) irradiance of each 3 RGB band',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY IRRADIANCE OF EACH 3 RGB BANDS JUST BELOW SEA ICE IN 2016
# --- umol photons m^-2 s^-1
ax=axs[4,0]
ax.plot(first_year,array1d_idoy_blue,
':',color='blue',lw=3,label='blue')
ax.plot(first_year,array1d_idoy_green,
'--',color='green',lw=3,label='green')
ax.plot(first_year,array1d_idoy_red,
'-.',color='red',lw=3,label='red')
ax.legend()
ax.set_ylabel('$E_\mathrm{i}$($\mathrm{0^-}$)(band)\n'
'($\mathrm{\mu}$mol photons $\mathrm{m^{-2}}$ $\mathrm{s^{-1}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,120)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(h) irradiance of each RGB band',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY IRRADIANCE OF EACH 3 RGB BANDS JUST BELOW OPEN WATER IN 2016
# --- umol photons m^-2 s^-1
ax=axs[4,1]
ax.plot(first_year,array1d_idoy_blueow,
':',color='blue',lw=3,label='blue')
ax.plot(first_year,array1d_idoy_greenow,
'--',color='green',lw=3,label='green')
ax.plot(first_year,array1d_idoy_redow,
'-.',color='red',lw=3,label='red')
ax.legend()
ax.set_ylabel('$E_\mathrm{w}$($\mathrm{0^-}$)(band)\n'
'($\mathrm{\mu}$mol photons $\mathrm{m^{-2}}$ $\mathrm{s^{-1}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,750)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(i) irradiance of each 3 RGB band',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY PAR JUST BELOW SEA ICE IN 2016
# --- umol photons m^-2 s^-1
ax=axs[5,0]
ax.plot(first_year,array1d_idoy_PARiceuEin,color='black',lw=3)
ax.set_ylabel('$I_\mathrm{i}$($\mathrm{0^-}$)\n'
'($\mathrm{\mu}$mol photons $\mathrm{m^{-2}}$ $\mathrm{s^{-1}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,120)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(j) PAR',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- DAILY PAR JUST BELOW OPEN WATER IN 2016
ax=axs[5,1]
ax.plot(first_year,array1d_idoy_PARowuEin,color='black',lw=3)
ax.set_ylabel('$I_\mathrm{w}$($\mathrm{0^-}$)\n'
'($\mathrm{\mu}$mol photons $\mathrm{m^{-2}}$ $\mathrm{s^{-1}}$)')
ax.set_xlim(first_year[0],first_year[-1])
ax.set_ylim(0,750)
ax.grid()
ax.set_xticks(locs)
ax.set_xticklabels(labels)
plt.text(0.50,1.10,
'(k) PAR',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes, fontsize=25)
# --- POSITION
axs[0,0].set_position( [0.10,0.75 ,0.89,0.22])
axs[0,1].set_position( [0.01,0.01 ,0.01,0.01])
axs[1,0].set_position( [0.10,0.58 ,0.40,0.10])
axs[1,1].set_position( [0.59,0.58 ,0.40,0.10])
axs[2,0].set_position( [0.10,0.44 ,0.40,0.10])
axs[2,1].set_position( [0.59,0.44 ,0.40,0.10])
axs[3,0].set_position( [0.10,0.30 ,0.40,0.10])
axs[3,1].set_position( [0.59,0.30 ,0.40,0.10])
axs[4,0].set_position( [0.10,0.16 ,0.40,0.10])
axs[4,1].set_position( [0.59,0.16 ,0.40,0.10])
axs[5,0].set_position( [0.10,0.02 ,0.40,0.10])
axs[5,1].set_position( [0.59,0.02 ,0.40,0.10])
with plt.style.context('mplstyles/par.mplstyle'):
# Plot
fig,axs=plt.subplots(nrows=6,ncols=2,figsize=(16, 16))
make_plots(axs)
# --- SAVE
plt.savefig('figures_supmat_progress/par.png',dpi=1000)
$m$ is the slope between the bloom initiation (I) and the bloom peak (P) in the log scale. A linear regression of the points between I and P was used.
m2010=blooms_df.loc[2010,['m']].item()
m2011=blooms_df.loc[2011,['m']].item()
m2015=blooms_df.loc[2015,['m']].item()
m2016=blooms_df.loc[2016,['m']].item()
blooms_df[['m']].describe()
| m | |
|---|---|
| count | 4.000000 |
| mean | 0.055439 |
| std | 0.017971 |
| min | 0.035786 |
| 25% | 0.043144 |
| 50% | 0.055373 |
| 75% | 0.067669 |
| max | 0.075224 |
def make_plots(ax):
box=dict(pad=5,alpha=0)
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
# --- RESOLUTE 2010
ax=axs[0]
h1=ax.plot(np.array(initiation2010_df.iT),
np.array(initiation2010_df.obsvintchl),
'o',color='red')
ax.set_yscale('log')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n(mg Chl $\mathrm{ m^{-2} }$)')
ax.set_ylim(0.5,1E3)
ax.grid()
# add slope between initiation and peak
h2=ax.plot(olsx2010,np.array(olsy2010),'-',color='red',lw=3)
h3=ax.plot(np.array(out2010_df.iT),
np.array(out2010_df.obsvintchl),
'o',color='black')
ax.text(90,110,'m={0:.3f}'.format(m2010),color='red',
fontsize=35)
plt.text(0.50,1.10,'(a) Resolute 2010',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- RESOLUTE 2011
ax=axs[1]
h1=ax.plot(np.array(initiation2011_df.iT),
np.array(initiation2011_df.obsvintchl),
'o',color='red')
ax.set_yscale('log')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n(mg Chl $\mathrm{ m^{-2} }$)')
ax.set_ylim(0.5,1E3)
ax.grid()
# add slope between initiation and peak
h2=ax.plot(olsx2011,np.array(olsy2011),'-',color='red',lw=3)
h3=ax.plot(np.array(out2011_df.iT),
np.array(out2011_df.obsvintchl),
'o',color='black')
ax.text(90,110,'m={0:.3f}'.format(m2011),color='red',
fontsize=35)
plt.text(0.50,1.10,'(b) Resolute 2011',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- QIKIQTARJUAQ 2015
ax=axs[2]
h1=ax.plot(np.array(initiation2015_df.iT),
np.array(initiation2015_df.obsvintchl),
'o',color='red')
ax.set_yscale('log')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n(mg Chl $\mathrm{ m^{-2} }$)')
ax.set_ylim(0.5,1E3)
ax.grid()
# add slope between initiation and peak
h2=ax.plot(olsx2015,np.array(olsy2015),'-',color='red',lw=3)
h3=ax.plot(np.array(out2015_df.iT),
np.array(out2015_df.obsvintchl),
'o',color='black')
ax.text(90,110,'m={0:.3f}'.format(m2015),color='red',
fontsize=35)
plt.text(0.50,1.10,'(c) Qikiqtarjuaq 2015',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- QIKIQTARJUAQ 2016
ax=axs[3]
h1=ax.plot(np.array(initiation2016_df.iT),
np.array(initiation2016_df.obsvintchl),
'o',color='red',
label='obs (0-100m) between initiation and peak')
ax.set_yscale('log')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylabel('Vertically integrated Chl $a$\n(mg Chl $\mathrm{ m^{-2} }$)')
ax.set_ylim(0.5,1E3)
ax.grid()
# add slope between initiation and peak
h2=ax.plot(olsx2016,np.array(olsy2016),'-',color='red',lw=3)
h3=ax.plot(np.array(out2016_df.iT),
np.array(out2016_df.obsvintchl),
'o',color='black',
label='obs (0-100m) before initiation and after peak')
ax.legend()
ax.text(90,110,'m={0:.3f}'.format(m2016),color='red',
fontsize=35)
plt.text(0.50,1.10,'(d) Qikiqtarjuaq 2016',
horizontalalignment = 'center', verticalalignment='center',
transform=ax.transAxes)
# --- ADDITIONAL X-AXIS
ax=axs[4]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.13,0.76 ,0.85,0.19])
axs[1].set_position( [0.13,0.52 ,0.85,0.19])
axs[2].set_position( [0.13,0.28 ,0.85,0.19])
axs[3].set_position( [0.13,0.04 ,0.85,0.19])
axs[4].set_position( [0.13,0.02 ,0.85,0.01])
with plt.style.context('mplstyles/m.mplstyle'):
# Plot
fig,axs=plt.subplots(nrows=5,ncols=1,figsize=(16.26,25.3))
make_plots(axs)
fig.align_ylabels(axs[:])
# --- SAVE
plt.savefig('figures_supmat_progress/m.png',dpi=1000)
Data is from CCGS Amundsen 2018, leg 2b. The location is shown on Fig. 1b (maps_bathymetry.png). Code by Raphaël Larouche and Maxime Benoît-Gagné.
def density_contour(ax, sgrid, tgrid):
"""
Plot of density contour.
:param ax:
:param sgrid:
:param tgrid:
:return:
"""
dens = gsw.sigma0(*np.meshgrid(sgrid, tgrid))
c_dens = ax.contour(sgrid, tgrid, dens, levels=np.arange(22, 29, 1), linestyles='dashed', colors='gray', linewidths=.8)
ax.clabel(c_dens, inline=True, fmt="%1.1f $\mathrm{kg\ m^{-3}}$", fontsize=10)
# Open data
path='data/DataS3_observations_CCGS_Amundsen/1802_ODV_truncated.csv'
data=pd.read_csv(path,header=[0],sep=',')
data=data.set_index("Station") # Set index of dataframe to station number
# Dictionary for station number
station_dic = {"station 1": 7,
"station 2": 9,
"station 3": 10,
"station 4": 12,
"station 5": 11,
"station 6": 13}
def make_plots(ax):
# Colormap for depth
colmap = mpl.cm.get_cmap("hot")
normalization = mpl.colors.Normalize(vmin=0.0, vmax=1100.0)
# Color each station
cl = {"station 1": "g",
"station 2": "b",
"station 3": "y",
"station 4": "c",
"station 5": "k",
"station 6": "m"}
# Linestyle each station
lstyle = {"station 1": "-",
"station 2": "-.",
"station 3": ":",
"station 4": "--",
"station 5": "-",
"station 6": "-."}
# Density lines
sgrid = np.linspace(29, 35, 100)
tgrid = np.linspace(-2, 2, 100)
density_contour(ax, sgrid, tgrid)
# Freeze curve
Sf = np.linspace(31, 34, 10)
Tf = gsw.CT_freezing(Sf, 0, 0)
ax.plot(Sf, Tf, "-", color="gray", linewidth=5)
# Loop
for key in station_dic.keys():
station_data=data.loc[station_dic[key]]
# Longitude, latitude
lon=station_data['Longitude [degrees_east]'].iloc[0]
lat=station_data['Latitude [degrees_north]'].iloc[0]
# Depth and normalization
depth=np.array(station_data["Depth [metres]"]).astype(list)
color_depth = [colmap(normalization(value)) for value in list(depth)]
# Salinity, conservative temperature data
salinity = np.array(station_data["Asal [g/kg]"]).astype(float)
cons_temp = np.array(station_data["Cont [deg C]"]).astype(float)
# Label
lab = "{0} ({1:.3f}$\degree$ W)".format(key, -1*lon)
# Plot data
if key == "station 5":
ax.plot(salinity, cons_temp, linewidth=5, color=cl[key], linestyle=lstyle[key], label=lab)
else:
ax.plot(salinity, cons_temp, linewidth=3, color=cl[key], linestyle=lstyle[key], label=lab)
# Figure TS diagram
ax.set_xlim(29, 35)
ax.set_ylim(-2, 2)
ax.set_yticks(np.arange(-1.5,2.5,0.5))
ax.set_ylabel('$\mathrm{ T_C\ (\degree C) }$')
ax.set_xlabel('$\mathrm{ S_A\ (g\ kg^{-1}) }$')
ax.legend(loc="best")
with plt.style.context('mplstyles/ts.mplstyle'):
# Plot
fig=plt.figure(figsize=(16, 16))
ax=fig.add_subplot(111)
make_plots(ax)
# --- SAVE
plt.savefig('figures_supmat_progress/ts.png',dpi=1000)
assert no3mod_df['doy'].equals(nutrients_df['doy']), \
"column doy in no3mod_df is different from column doy in nutrients_df"
assert no3mod_df['depth'].equals(nutrients_df['depth']), \
"column depth in no3mod_df is different from column depth in nutrients_df"
s=nutrients_df['no3'].copy()
s[s==0]=np.NaN
array1d_iobs_logno3obs=np.log10(s)
array1d_iobs_logno3mod =np.log10(no3mod_df['no3mod'])
refnstd=np.nanstd(array1d_iobs_logno3obs,ddof=1) \
/np.nanstd(array1d_iobs_logno3obs,ddof=1)
nstdlogno3=np.nanstd(array1d_iobs_logno3mod,ddof=1) \
/np.nanstd(array1d_iobs_logno3obs,ddof=1)
corrcoeflogno3=(pd.DataFrame(
{'array1d_iobs_logno3obs':array1d_iobs_logno3obs,
'array1d_iobs_logno3mod':array1d_iobs_logno3mod}) \
.corr()).iloc[0,1]
assert sioh4mod_df['doy'].equals(nutrients_df['doy']), \
"column doy in sioh4mod_df is different from column doy in nutrients_df"
assert sioh4mod_df['depth'].equals(nutrients_df['depth']), \
"column depth in sioh4mod_df is different from column depth in nutrients_df"
s=nutrients_df['sioh4'].copy()
s[s==0]=np.NaN
array1d_iobs_logsioh4obs=np.log10(s)
array1d_iobs_logsioh4mod =np.log10(sioh4mod_df['sioh4mod'])
refnstd=np.nanstd(array1d_iobs_logsioh4obs,ddof=1) \
/np.nanstd(array1d_iobs_logsioh4obs,ddof=1)
nstdlogsioh4=np.nanstd(array1d_iobs_logsioh4mod,ddof=1) \
/np.nanstd(array1d_iobs_logsioh4obs,ddof=1)
corrcoeflogsioh4=(pd.DataFrame(
{'array1d_iobs_logsioh4obs':array1d_iobs_logsioh4obs,
'array1d_iobs_logsioh4mod':array1d_iobs_logsioh4mod}) \
.corr()).iloc[0,1]
assert po4mod_df['doy'].equals(nutrients_df['doy']), \
"column doy in po4mod_df is different from column doy in nutrients_df"
assert po4mod_df['depth'].equals(nutrients_df['depth']), \
"column depth in po4mod_df is different from column depth in nutrients_df"
s=nutrients_df['po4'].copy()
s[s==0]=np.NaN
array1d_iobs_logpo4obs=np.log10(s)
array1d_iobs_logpo4mod =np.log10(po4mod_df['po4mod'])
refnstd=np.nanstd(array1d_iobs_logpo4obs,ddof=1) \
/np.nanstd(array1d_iobs_logpo4obs,ddof=1)
nstdlogpo4=np.nanstd(array1d_iobs_logpo4mod,ddof=1) \
/np.nanstd(array1d_iobs_logpo4obs,ddof=1)
corrcoeflogpo4=(pd.DataFrame(
{'array1d_iobs_logpo4obs':array1d_iobs_logpo4obs,
'array1d_iobs_logpo4mod':array1d_iobs_logpo4mod}) \
.corr()).iloc[0,1]
array1d_iobs_logchlHPLC=np.log10(chlHPLC_df['chlHPLC'])
array1d_iobs_logchlmod =np.log10(chlmod_df['chlmod'])
refnstd=np.std(array1d_iobs_logchlHPLC,ddof=1) \
/np.std(array1d_iobs_logchlHPLC,ddof=1)
nstdlogchl=array1d_iobs_logchlmod.std(ddof=1) \
/array1d_iobs_logchlHPLC.std(ddof=1)
corrcoeflogchl=(pd.DataFrame(
{'array1d_iobs_logchlHPLC':array1d_iobs_logchlHPLC,
'array1d_iobs_logchlmod':array1d_iobs_logchlmod}) \
.corr()).iloc[0,1]
def make_plots(axs):
plt.axis('off')
dia=TaylorDiagram.TaylorDiagram(refnstd,
fig=fig,
rect=111,
label='Reference',
srange=(0,2.5),
extend=False)
dia.add_sample(nstdlogchl,
corrcoeflogchl,
marker='o',
ms=10,
ls='',
mfc='red',
mec='red')
dia.add_sample(nstdlogno3,
corrcoeflogno3,
marker='o',
ms=10,
ls='',
mfc='red',
mec='red')
dia.add_sample(nstdlogsioh4,
corrcoeflogsioh4,
marker='o',
ms=10,
ls='',
mfc='red',
mec='red')
dia.add_sample(nstdlogpo4,
corrcoeflogpo4,
marker='o',
ms=10,
ls='',
mfc='red',
mec='red')
dia.add_grid()
contours=dia.add_contours(colors='0.5')
plt.clabel(contours, inline=1, fontsize=10, fmt='%.2f')
plt.text(0.39,0.24,'CHL',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, color='red', fontsize=20)
plt.text(0.20,0.11,'NO3',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, color='red', fontsize=20)
plt.text(0.33,0.16,'SIL',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, color='red', fontsize=20)
plt.text(0.27,0.18,'PO4',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, color='red', fontsize=20)
plt.text(0.41,0.03,'REF',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, color='red', fontsize=20)
plt.tight_layout()
with plt.style.context('mplstyles/taylor.mplstyle'):
# Plot
fig=plt.figure(figsize=(15,10))
ax=fig.add_subplot(111)
make_plots(axs)
# --- SAVE
plt.savefig('figures_supmat_progress/taylor.png',dpi=1000)
def make_plots(axs):
box=dict(pad=5,alpha=0)
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
# --- ICE CONCENTRATION
ax=axs[0]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)',bbox=box)
ax.set_ylim(0,100)
ax.grid()
ax.legend(loc='upper left',bbox_to_anchor=(0.05,1.00))
plt.text(0.03,0.6,'(a)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- VERTICALLY INTEGRATED 0-40m
ax=axs[1]
h1=ax.plot(first_year,array1d_iT_modvint40mchl,
'-',lw=3,color='black',label='model Chl')
h2=ax.plot(first_year,array1d_iT_obsvint40mchlHPLC,
'o',markersize=6,color='black',label='obs Chl')
ax.set_yscale('log')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylim(2E-1,3E2)
ax.set_ylabel('Vertically integrated Chl $a$\n'
'(0-40m, $\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.grid()
axt = ax.twinx()
# g C m^-2 -> mg C m^-2
h3=axt.plot(first_year,array1d_iT_modvint40mphytobio*1E3,
'--',lw=3,color='red',label='model C')
axt.set_yscale('log')
axt.set_ylim(*[100*l for l in ax.get_ylim()])
axt.set_ylabel('Vertically integrated\n'
'phyto biomass\n'
'(0-40m, $\mathrm{ mg\ C\ m^{-2} }$)')
# ref: https://stackoverflow.com/questions/5484922/secondary-axis-with-twinx-how-to-add-to-legend
# answer by Paul
# lns=h1+h2+h3
lns=h3+h1+h2
labs = [l.get_label() for l in lns]
ax.legend(lns,labs,loc='upper left',bbox_to_anchor=(0.05,1.00))
plt.text(0.03,0.93,'(b)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- VERTICALLY INTEGRATED 0-100m
ax=axs[2]
h1=ax.plot(first_year,array1d_iT_modvintchl,
'-',lw=3,color='black',label='model Chl')
h2=ax.plot(first_year,array1d_iT_obsvintchlHPLC,
'o',markersize=6,color='black',label='obs Chl')
ax.set_yscale('log')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylim(2E-1,3E2)
ax.set_ylabel('Vertically integrated Chl $a$\n'
'(0-100m, $\mathrm{ mg\ Chl\ m^{-2} }$)')
ax.grid()
axt = ax.twinx()
# g C m^-2 -> mg C m^-2
h3=axt.plot(first_year,array1d_iT_modvintphytobio*1E3,
'--',lw=3,color='red',label='model C')
axt.set_yscale('log')
axt.set_ylim(*[100*l for l in ax.get_ylim()])
axt.set_ylabel('Vertically integrated\n'
'phyto biomass\n'
'(0-100m, $\mathrm{ mg\ C\ m^{-2} }$)')
plt.text(0.03,0.93,'(c)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- MEAN 0-40m
ax=axs[3]
h1=ax.plot(first_year,array1d_iT_modmeanchl,
'-',lw=3,color='black',label='model Chl')
h2=ax.plot(first_year,array1d_iT_obsmeanchlHPLC,
'o',markersize=6,color='black',label='obs Chl')
ax.set_yscale('log')
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_ylim(7E-3,5E0)
ax.set_ylabel('Mean Chl $a$\n'
'(0-40m, $\mathrm{ mg\ Chl\ m^{-3} }$)')
ax.grid()
axt = ax.twinx()
# g C m^-3 -> mg C m^-3
h3=axt.plot(first_year,array1d_iT_modmeanphytobio*1E3,
'--',lw=3,color='red',label='model C')
axt.set_yscale('log')
axt.set_ylim(*[100*l for l in ax.get_ylim()])
axt.set_ylabel('Mean\n'
'phyto biomass\n'
'(0-40m, $\mathrm{ mg\ C\ m^{-3} }$)')
plt.text(0.03,0.93,'(d)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[4]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.13,0.92 ,0.74,0.04])
axs[1].set_position( [0.13,0.63 ,0.74,0.27])
axs[2].set_position( [0.13,0.34 ,0.74,0.27])
axs[3].set_position( [0.13,0.05 ,0.74,0.27])
axs[4].set_position( [0.13,0.03 ,0.74,0.01])
with plt.style.context('mplstyles/vint40m.vint100m.mean40m.mplstyle'):
# Plot
fig,axs=plt.subplots(nrows=5,ncols=1,figsize=(18.4,17.8))
make_plots(axs)
fig.align_ylabels(axs[:])
# --- SAVE
plt.savefig('figures_supmat_progress/chl.bio.png',dpi=1000)
def make_plots(axs):
locs=np.array([0, 31, 59, 90, 120, 151,
181, 212, 243, 273, 304, 334])
labels=('Jan-1','Feb-1','Mar-1','Apr-1','May-1','Jun-1',
'Jul-1','Aug-1','Sep-1','Oct-1','Nov-1','Dec-1')
# --- ICE CONCENTRATIONS
ax=axs[0]
h1=ax.plot(first_year,ice*100,'k-',lw=2.5,label='model')
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks([0,50,100])
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Ice conc.\n (%)')
mpl.axis.Axis.set_label_coords(ax.yaxis,-0.06,0.20)
ax.set_ylim(0,100)
ax.grid()
plt.text(0.03,0.60,'(a)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- ACCUMULATION RATE
ax=axs[1]
h1=ax.plot(first_year,array1d_iT_r2016,'-',color='blue',lw=3)
ax.set_xlim(first_year[0],first_year[-1])
xlims=ax.get_xlim()
ax.set_xticks(locs)
ax.set_xticklabels([])
ax.set_yticks(np.arange(-0.5,1.0,0.5))
ax.yaxis.set_ticks_position('both')
ax.set_ylim(-0.5,0.5)
ax.set_ylabel('Accumulation rate\n($\mathrm{ d^{-1} }$)')
ax.plot(xlims,[0,0],'k--',color='gray')
ax.grid()
plt.text(0.03,0.93,'(b)',
horizontalalignment = 'center',verticalalignment='center',
transform=ax.transAxes, fontsize=30)
# adding the vertical line at start of bloom development
ax.axvline(tI-1,color='green',linestyle='--')
# adding the vertical line at 15% of observed ice
ax.axvline(obs_breakup,color='k',linestyle='--')
# --- ADDITIONAL X-AXIS
ax=axs[2]
ax.set_xlim(xlims)
ax.set_xticks(locs)
ax.set_xticklabels(labels)
ax.axes.get_yaxis().set_visible(False)
ax.spines['top'].set_visible(False)
# --- POSITION
axs[0].set_position( [0.13,0.88 ,0.85,0.10])
axs[1].set_position( [0.13,0.10 ,0.85,0.73])
axs[2].set_position( [0.13,0.07 ,0.85,0.01])
with plt.style.context('mplstyles/exp.mplstyle'):
# Plot
fig,axs=plt.subplots(3,1,figsize=(16,6.6))
make_plots(axs)
fig.align_ylabels(axs[:])
# --- SAVE
plt.savefig('figures_supmat_progress/r.png',dpi=1000)